Custom types for decompositions

[Andlon]

Rulinalg has implementations of a number of decompositions. Currently these are returned as Result<tuple, Error>, where tuple is a tuple of Matrix<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 of b. 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):

let qr = A.qr();
for b in right_hand_sides {
    let x = qr.solve(b);
    do_something(x);
}

This is possible if qr() returns a struct QR that has a specialized solve method, as well as accessors to the decomposed parts Q and R. Similarly for LU. 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!

[AtheMathmo]

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.

[Andlon]

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.

[Andlon]

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.

[AtheMathmo]

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.

[Andlon]

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

[AtheMathmo]

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.

[Andlon]

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

[AtheMathmo]

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!

[AtheMathmo]

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.

[Andlon]

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.

[AtheMathmo]

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.

[Andlon]

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?

[AtheMathmo]

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

[Andlon]

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.