-
Notifications
You must be signed in to change notification settings - Fork 60
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Custom types for decompositions #40
Comments
Thanks for writing this up -I think it's a great idea! It certainly makes the decompositions easier to use in some cases and makes the signatures of some functions nicer. I worry that the api could be a little too verbose though, for example I could see things like the following being inevitable: let (q, r) = A.qr().unwrap().get_parts(); I still think that overall it is a lot cleaner - but documentation will be essential to keep things obvious to the end user. I also wonder if there is a way to use this to help us access the different types of algorithms for each decomposition. For example for SVD we may only want the singular values which takes less time. Or even eigenvalues/vector algorithms. Right now the only way we can do this is by providing different functions - but perhaps there is a tidier way? I'm not sure I'll have time to work on this for a while but I'd be really happy if someone else had the time to pick it up. |
Yeah, the slightly more verbose syntax if you just want the Q and R matrices is unavoidable, but is a worthwhile trade-off in my opinion. I've been thinking about your SVD example today. I'm not sure if I've interpreted you correctly, but let's see where it takes us. The way I see it, there are two ways forward:
Let's toy with the idea of decompositions being expressions. Using your example, say we first are only interested in the singular values. Then we could do something like
What do we do here? Well, I guess first we computed the singular values in the most efficient way possible, but then we're asked to compute the pseudo-inverse, so the data we already have is insufficient, and we'd have to re-compute the SVD. I once had a professor who referred to the SVD as the "swiss-army knife of linear algebra". With this comes the implication that the SVD is not likely the most efficient approach for any one of its many applications, but its flexible and robust, and in many cases still offers decent performance. From a user's perspective, it's a quick and convenient way to obtain the right results. If you really need to obtain those results in a faster way, you can always resort to very specialized algorithms. My gut feeling tells me that the right approach here is to keep it simple. For the user, this means no surprises. The user asks us to compute the SVD, so we do it right then and there. Once it's obtained, he can freely use the results of the full SVD. If he only needs, say, the largest singular value or so, there are better suited algorithms, but even then it might be useful to the user for a first implementation (before he starts caring about performance and scale). There may be more clever approaches that I haven't thought of. I'm also not too well versed in Rust yet. |
Btw, if we define
but that assumes that |
I actually really like this lazy evaluation approach. It feels very rust-like and has the advantage that it lets us compute things in the most efficient manor. The code is a little more verbose but the API will be easier to navigate I think - the user doesn't have to find the I agree that this does add some complication - especially for users who aren't so familiar with Rust. But I think that adding a I'd definitely likely to at least explore this (perhaps for eigen decomp first) to see how it feels. |
Could you perhaps sketch out some pseudocode that demonstrates roughly what you have in mind? I'm not sure if I'm quite following! And yeah, I'm all for exploring :) |
Sure! Let me know if things still aren't clear (or if you see any issues or just plain disagree): use libnum::Zero;
/// Dummy Eigen decomp
#[must_use = "Eigen decomp is lazy and must be used"]
pub struct EigenDecomp<T> {
// Stuff would go here
}
impl<T: Clone + Zero> EigenDecomp<T> {
/// Consumes decomp to get eigen values
pub fn values(self) -> Matrix<T> {
// Compute the eigen values
}
/// Consumes decomp to get eigen vectors
pub fn vectors(self) -> Matrix<T> {
// Compute the eigen vectors
}
}
impl<T: Zero> Matrix<T> {
/// Returns a lazy EigenDecomp computer
pub fn eigdecomp(&self) -> EigenDecomp<T> {
EigenDecomp {
// Stuff would go here (probably take ownership of self..)
}
}
} The above is a dummy implementation of eigen decomposition. In practice there would probably be some differences in the api too. The above would mean that the following code would produce a compiler warning as the return value is not used (similar to let eigs = A.eigdecomp(); And to get the eigen values the user would do this: let values = A.eigdecomp().values(); I think the advantage is that there is one entry point for computing eigen-things whereas now we have
And other variants exist too. |
Thanks for elaborating! I think I see what you mean now. If I understand you correctly, the key point of your proposal is that we get to collect related functionality within thinner "expression" objects, from which one can further perform the relevant computations. This certainly has merits, I agree! One problem with your example of eigen decomposition, however, is that you can compute eigenvalues for a matrix for which an eigendecomposition does not exist. Hence, if you seek to compute eigenvalues for such a matrix, it would be odd/confusing to go through the eigendecomposition in the API. However, this is easily remedied:
Perhaps Similarly, for SVD:
What about QR/LU(P)/Cholesky? Perhaps here we can go straight to the decomposition, since there's no obvious need (that I can think of) for not doing it. E.g.
I think the above looks pretty good and intuitive. It essentially combines our two versions of the proposal, hopefully taking the best of both. I'm fairly tired right now, so I might have missed something. By the way, as for the ordering of the singular values... It's a bit comparing apples to oranges, since one is measured in flops and the other is essentially memory accesses, but SVD for a square matrix is approx. I actually wanted to try to implement the sorting of the singular values, as I saw the issue was up for grabs and I figured it was a relatively small thing that I could use for a first contribution (assuming std sort). Maybe next week. Or maybe not :) |
Yeah I think this proposal makes a good amount of sense! There is a little disparity as you say with QR/Cholesky/etc. but I don't see this being a significant obstacle. One last point I'd make is that to me this makes it a lot easier to provide bindings to LAPACK for these decomposition algorithms. Right now I haven't been able to find a good way to do so but using new modules the thinner "expression" objects will make it easy to implement the various LAPACK algorithms under the hood and expose them cleanly. At least to me this seems like a step in the right direction.
For your SVD comments I'll address them in issue #17 . I think there are some valuable points that will be useful for whoever implements this.
Hopefully you! |
I'd also like to point out that both the Q and R matrices of QR can be stored in a single matrix (see for example here). This is also the case for LU decomposition. A benefit of encapsulating the decompositions in custom types is therefore that we can consume One might think that first storing it this way and then reconstructing them might be a performance problem if the user only wants |
I agree that this is definitely worth exploring! Will have more to say on this later but I think for now we should tackle eigen and svd as they are slightly simpler to convert. |
I think we could actually reimplement all decompositions as structs now, but keep the current functionality. I.e. for QR, we internally just store the current What I'm suggesting is to overhaul the decompositions API in the following way:
Number 1 can be done simultaneously for all decompositions in its own PR. This will make it easier to apply subsequent changes, and people can more easily work on individual decompositions at the same time. I can volounteer to try to do this fairly soon (this week, hopefully?). Number 2 can then perhaps be done on a per-decomposition basis, and number 3 in the longer run, depending on our priorities. One problem with this approach, is that for some decompositions, there is not so much utility in having a specialized struct. For example, I don't know of any special applications of the upper/lower Hessenberg decomposition, and so having a specialized struct here is quite redundant (unless such an application exists). We could keep these decompositions as they are, but it feels a bit inconsistent. However, that might not be a problem. What do you think? |
This sounds like a good idea. In response to your points:
I agree that tackling things like upper/lower Hessenberg are a little more overkill - for now we can keep those utility functions (givens rotations, etc.) in the decomposition root module. If you do tackle this feel free to rewrite totally or work from #44 |
I think perhaps we should think a little independently about Eigen/SVD and LU/QR. In the first case, it makes sense as you say to have a "thin" entry point, although in effect it is really just a grouping of functionality (it doesn't actually give us any new functionality over just having these functions directly in In the second case, I'm suggesting that For SVD, you might even have multiple "full" decompositions, so it would make sense to be able to call i.e. So, to summarize, |
@AtheMathmo: I can see if I have some time this week and look into creating a preliminary PR. Might make it easier for us if we have something more tangible to work with. |
Just commenting to say that I think the above makes a lot of sense and is also what I was thinking. I'll be fairly unavailable next week (on vacation) but I look forward to seeing the PR when you (and I) do get time! |
Great, just wanted to clarify what I had in mind. Thanks for the headsup, and enjoy your vacation! :) |
Having considered this a bit longer, I'm actually not sold on the utility of the Even the argument for LAPACK does not really hold, as if I understand correctly, it's just a manner of organizing the implementations in different modules and using |
It may be that I cannot convince you - but I do think there is still some value in this structure. I agree that certainly we can achieve much of the same by separating into modules and using cfg directives. But my argument involving LAPACK was more to help collect the many different lapack routines that can be used for computing eigen values. I was trying to argue that we can put the entry point for all of these functions within the I'm happy to see a pull request removing the |
So, if I understand you correctly, then the idea is to expose direct wrappers around the Lapack functions, such as for example In that case, I can understand the motivation much better. I thought you wanted to expose the same API, but with different backends (i.e. our own vs LAPACK), but you actually want to expose an extended API covering additional LAPACK methods if the optional LAPACK support for EDIT: I just noticed a lot of the functions are the same, but for different types, so my particular example of |
I was initially thinking that we would directly expose some routines and then have higher level functions for more common ones. So But honestly, introducing LAPACK is a pretty huge step and I expect our approach to change as we move forwards. I have no experience using LAPACK directly and deciding how the user should interact with it seems tricky! |
I'm currently more of the following opinion, but I can quite possibly be convinced: If users really know what they are doing and want direct LAPACK access, they would probably be better off with using LAPACK directly or a dedicated crate whose purpose is essentially wrapping LAPACK. I think perhaps However, in order to get to this (completely arbitrary and made-up) 99%, we'd have to provide access to some common specializations. We could maybe expose some specialized functionality using more idiomatic rust. To use your example, consider tridiagonal symmetric matrices. In this case, we could imagine something like (note that I'm just throwing out ideas here, not suggesting this is the best way to solve it)
with usage:
You'll notice this is similar to the
|
Anyway, I wanted to start on reorganizing all the decompositions (i.e. in terms of splitting into multiple files, no change in API or functionality), but I was unsure whether to branch off |
I still need to read through your above comment properly (but am currently agreeing with it mostly). As for the splitting into multiple files - I think doing this in a fresh PR is the better approach. Thanks for all your time put into this - I really appreciate the thought and detailed commentary you've provided! |
What is the status of this? It seems like what is suggested has been implemented in Edit: Actually it is implemented on master, just not yet in the version on crates.io. |
@vks: You are right, there's a lot of stuff in the pipeline which is not yet on crates.io. In particular, on master there's also |
Rulinalg has implementations of a number of decompositions. Currently these are returned as
Result<tuple, Error>
, wheretuple
is a tuple ofMatrix<T>
. This is in itself useful, as it gives the user easy access to the decomposed parts. However, often you are not necessarily so interested in the decomposed parts themselves, but rather what you can do with them.For example, in the case of LU or QR and a linear system
Ax = b
, you may want to perform the decomposition once, and then use it many times over to solve for different instances ofb
. In this case, as a user, you don't want to have to perform the transposition of Q and the solution of the triangular system. Rather, it would be convenient if you could do something along the lines of (ignoring error handling for to make my point more clear):This is possible if
qr()
returns a structQR
that has a specializedsolve
method, as well as accessors to the decomposed partsQ
andR
. Similarly forLU
. In my opinion, you would not lose anything over the current API, but you'd make it easier to use, and it could also open up for optimizations in how the decompositions are stored, and how they are used to perform certain actions.For SVD, you could have convenient methods to retrieve the singular values, perhaps for least squares solution etc. In addition to accessors to the decomposed parts, of course.
Eigen (C++ library) does something like this. I'll be quick to add that I do not want
rulinalg
to mirror Eigen's API too closely, as it is... perhaps overly complex at times.Do you think this style of API could be suitable for
rulinalg
? I'd be interested to hear your opinions!The text was updated successfully, but these errors were encountered: