Custom types for decompositions about rulinalg HOT 26 OPEN

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.

from rulinalg.

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:

• A decomposition is a simple struct containing the necessary data, with some extra convenient functions for typical applications of that decomposition, such as solving a linear system repeatedly much more efficiently.
• A decomposition is some kind of "expression", meaning it is not evaluated directly, but rather lazily in some way, depending on the next actions of the user.

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 A.svd().singular_values(). Great, because the decomposition is lazy, it was able to compute the singular values in the most efficient way possible. Now, say we want to obtain the singular values and compute the pseudo-inverse of A. We could do:

 // Note that we need the mut if we want to use the lazy approach, which is unfortunate
 let mut svd = A.svd().unwrap();
 let s = svd.singular_values();
 let p = svd.pseudo_inverse();

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.

from rulinalg.

Btw, if we define QR as a tuple struct, we should be able to do:

let QR(q, r) = A.qr().unwrap();

but that assumes that q and r correspond to public fields in the tuple struct. We really don't want to make them public as then we cannot ensure the necessary invariants, since the user is able to change them. I don't know Rust well enough to say if it's possible to customize the destructuring process in such a way that we can let the fields remain private? Not to my knowledge.

from rulinalg.

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 eigenvalues function to know it is better to use that than eigendecomp if they only need the values.

I agree that this does add some complication - especially for users who aren't so familiar with Rust. But I think that adding a must_use attribute would help with ensuring that the object is consumed.

I'd definitely likely to at least explore this (perhaps for eigen decomp first) to see how it feels.

from rulinalg.

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 :)

from rulinalg.

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 Result).

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 eigenvalues() and eigendecomp(). For SVD this is an even bigger problem as we could have things like:

• Just singular values
• Just singular values un-ordered (If this is the same complexity I'd love to know as I don't have it implemented this way right now...)
• Full decomposition
• Full decomposition with zero blocks

And other variants exist too.

from rulinalg.

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:

let values = A.eigen().values().unwrap();
let vectors = A.eigen().vectors().unwrap();
let decomposition = A.eigen().decomposition().unwrap();
let values_from_decomp = decomposition.values();
let vectors_from_decomp = decomposition.vectors();

Perhaps eigen().vectors() is unnecessary, since - at least as far as I know - you need to compute the eigenvalues before you compute the eigenvectors anyway, and as such you're probably better off making the full decomposition, since then you can also extract the values. Although the term decomposition is only valid if it is diagonalizable. Eh, let's eschew this detail for now.

Similarly, for SVD:

let values = A.svd().values().unwrap();

// Perform thin, full, compact, truncated decompositions
let thin = A.svd().thin().unwrap();
let full = A.svd().full().unwrap();
// ...
let values_from_full = full.values();
let left_singular_vectors = full.left();

// Immutably borrow the matrices
let (u, sigma, v) = full.parts(); // Not sure if 'parts' is the best word

// Take ownership of the matrices and consume `full`
let (u, sigma, v) = full.take_parts(); 

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.

let qr = A.qr().unwrap();
// Can solve for x in I think `O(n^2)` now iirc
let x = qr.solve(b); // If x is square
let y = qr.least_squares(b); // If m > n
let (q, r) = qr.parts(); // Borrow
let (q, r) = qr.take_parts(); // Take ownership

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. O(n^3) iirc. You can sort the eigenvalues in O(n log n) time. If using the eigenvalues as the sorting key, you can simultaneously get the sort order of the eigenvectors. If we allocate a new matrix, you can simply loop over all matrix elements, using the obtained sort order of the vectors to obtain each element. Hence, this is merely O(n^2). So, in this case, asymptotically, the sort should be relatively negligible. In practice, sorting is also very fast (or can be very fast, I'm not familiar with Rust's sort in particular), so I expect it also to be more or less negligible. Since you know the final position of each vector, I believe you can also perform the operation in-place in the same amount of time by swapping columns, though I'm too tired to verify. By the way, for this case, it might be a bit of a disadvantage that rulinalg uses row-major matrices.

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 :)

from rulinalg.

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.

// This will use rulinalg's eigvalue algorithm by default
// Or LAPACK's if the lapack feature is enabled
let values = A.eigen().values().unwrap();

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.

I actually wanted to try to implement the sorting of the singular values

Hopefully you!

from rulinalg.

I created #44 to explore what this looks like for eigen methods. If you have time to review this @Andlon I'd appreciate it. I will likely try to do the same for SVD next and maybe explore exposing lapack here too.

from rulinalg.

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 self on a matrix and construct a QR/LU decomposition using no more storage (except temporary variables in the algorithm itself) than the original matrix. We can expose an interface to let the user recover the separate Q and R matrices if he desires.

One might think that first storing it this way and then reconstructing them might be a performance problem if the user only wants Q and R directly, but I think it would not be an issue. Recall that both QR and LU for a square matrix are O(n^3), while the reconstruction is O(n^2) (I think?). In the case of LU, this is very fast, as it merely involves copying half the elements of the matrix and zeroing the rest. In the QR case I'm not sure exactly how fast the reconstruction is (recovering Q is a bit more work), but I think it should still be more or less negligible in comparison to the actual decomposition. Of course, we'd need to bench this.

from rulinalg.

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.

from rulinalg.

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 Q and R matrices that are currently generated. That way we can make relevant optimizations (storage etc.) without having to change the API (ideally). Does that make sense?

What I'm suggesting is to overhaul the decompositions API in the following way:

1. Split each decomposition into its own internal module/file (but don't expose this to the end-user)
2. Rewrite each decomposition to return a struct that encapsulates the parts of the decomposition, with an API to access the individual parts.
3. Add specialized features and improve optimality of storage to these decomposition structs. This would include adding specialized methods like solve() or least_squares(), depending on the decomposition in question.

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?

from rulinalg.

This sounds like a good idea. In response to your points:

1. We should definitely do this - however for documentation we will need to expose these to the end-user. I think...
2. I also agree with this - but I'm not sure how consistent we can have the API across each decomposition. For example part of the benefit of this approach, I believe, is to let the user choose optimal algorithms within a thinner entry point. e.g. for eigen decomposition they can compute eigen values and/or vectors. I'd rather not fill the struct with Option fields to handle this. Does that make sense?
3. Totally agree with this - and that it should come later.

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

from rulinalg.

I also agree with this - but I'm not sure how consistent we can have the API across each decomposition. For example part of the benefit of this approach, I believe, is to let the user choose optimal algorithms within a thinner entry point. e.g. for eigen decomposition they can compute eigen values and/or vectors. I'd rather not fill the struct with Option fields to handle this. Does that make sense?

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 Matrix, but hopefully it makes for a better user experience).

In the second case, I'm suggesting that matrix.qr() would perform the QR decomposition immediately, returning a QR struct that contains the result of the decomposition. This struct comes with implementations of common applications of the decomposition, such as i.e. lu().solve(b) or qr().least_squares(b). This lets the user easily reuse the decomposition for i.e. solving with multiple right-hand sides, and everything can be taken care of by rulinalg in an optimal manner.

For SVD, you might even have multiple "full" decompositions, so it would make sense to be able to call i.e. mat.svd().full() or mat.svd().thin() or various functions for low-rank approximations etc. The result of each of these could be similar to the QR struct, as it holds the results of the decomposition, and you can use it for various applications (i.e. with SVD you can also have least_squares()). The SVD has a lot of applications, and this could be a natural place to add such functionality.

So, to summarize, matrix.svd() or matrix.eigen() would only give a thin wrapper, whereas matrix.qr() or matrix.lu() would give "populated" structs. Further methods on eigen or svd could give "full structs", like those of QR or LU.

from rulinalg.

@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.

from rulinalg.

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!

from rulinalg.

Great, just wanted to clarify what I had in mind. Thanks for the headsup, and enjoy your vacation! :)

from rulinalg.

Having considered this a bit longer, I'm actually not sold on the utility of the eigen() organization. I really cannot see any clear benefits to it compared to simply matrix.eigen() to compute an eigendecomposition and matrix.eigenvalues() to compute only the eigenvalues.

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 cfg directives to switch between them.

from rulinalg.

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 Eigen struct. Hoping that this would make things a little easier for the user - of course using a module which we expose to the user may achieve the same.

I'm happy to see a pull request removing the eigen() entry point - maybe it is a little convoluted :)

from rulinalg.

So, if I understand you correctly, then the idea is to expose direct wrappers around the Lapack functions, such as for example mat.eigen().sorgtr()? (Just to take an arbtirary 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 rulinalg has been enabled...?

EDIT: I just noticed a lot of the functions are the same, but for different types, so my particular example of sorgtr() was not the best, as it's specifically for f32. We'd of course like to use the type system for this .

from rulinalg.

I was initially thinking that we would directly expose some routines and then have higher level functions for more common ones. So mat.eigen().pteqr() would be exposed so users can have a more efficient call for tridiagonal symmetric matrices. And we would want to also have mat.eigen().values() use LAPACK routines if the feature flag is set.

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!

from rulinalg.

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 rulinalg could be a performant, easy-to-use alternative which would be enough for 99% of users. Only the most extreme power users would likely need direct LAPACK access.

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)

pub struct SymmetricTridiagonal<M>(M);

impl<T> SymmetricTridiagonal<Matrix<T>> {
    // Under the hood we can easily use LAPLACK or whatever we want
    pub fn eigenvalues(&self) -> Vector<T> { ... }
    pub fn solve(&self, b: Vector<T>) -> Vector<T> { ... }

    // Whatever else you may want to do with a symmetric tridiagonal matrix
}

with usage:

// from_matrix checks that the invariant (symmetric + tridiagonal) holds,
// and returns Option or Result
let mat = SymmetricTridiagonal::from_matrix(other_mat).unwrap();

// from_matrix_unchecked does no checking at all. Maybe it should even be unsafe
// to make sure our users think about this properly? I.e. it should only be used
// after having made sure the safe alternative is not a bottleneck
let mat = SymmetricTridiagonal::from_matrix_unchecked(other_mat);

let values = mat.eigenvalues();

You'll notice this is similar to the LowerTriangular suggestion I made earlier in #41. There may be better ways, but these wrappers are the best I've come up with so far... Especially because this also solves the problem of "approximately fulfilled" invariants, for example if the matrix is only symmetric tridiagonal up to some low tolerance due to numerical noise, you could do

// Consider any value with absolute value smaller than 1e-8 to be zero
let mat = SymmetricTridiagonal::from_matrix_with_tol(other_mat, 1e-8).unwrap();

from rulinalg.

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 master, or go from your eigen() PR. What's your opinion?

from rulinalg.

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!

from rulinalg.

What is the status of this? It seems like what is suggested has been implemented in PartialPivLu, but not for other decompositions, like Cholesky.

Edit: Actually it is implemented on master, just not yet in the version on crates.io.

from rulinalg.

@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 FullPivLu, Cholesky, HouseholderQr and in a not-yet-merged PR (#179) a rewrite of the Hessenberg decomposition. I've been very busy lately, but I hope to pick up the pace and in the not-too-distant future provide a real schur decomposition (which should allow us to provide more robust eigenvalue/eigenvector computation) and a reworked SVD.

from rulinalg.