Polar decomposition

[LucasCampos]

Hello. As detailed in this thread in the user forum, I am implementing the matricial polar decomposition in nalgebra.

The bulk of the code is ready. However, I am having some issues with the traits needed to properly compile the code.

Cleanup and further implementations are still in the 'to do' list.

Any pointers on the specific type traits needed would be appreciated.

The compile errors are

error[E0277]: cannot multiply `base::matrix::Matrix<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic>>` to `base::matrix::Matrix<N, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<N, base::dimension::Dynamic, base::dimension::Dynamic>>`
  --> src/linalg/polar.rs:76:39
   |
76 |                     Some(v_t.adjoint()*DMatrix::from_diagonal(&svd.singular_values)*v_t)
   |                                       ^ no implementation for `base::matrix::Matrix<N, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<N, base::dimension::Dynamic, base::dimension::Dynamic>> * base::matrix::Matrix<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic>>`
   |
   = help: the trait `std::ops::Mul<base::matrix::Matrix<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic>>>` is not implemented for `base::matrix::Matrix<N, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<N, base::dimension::Dynamic, base::dimension::Dynamic>>`
   = help: consider adding a `where base::matrix::Matrix<N, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<N, base::dimension::Dynamic, base::dimension::Dynamic>>: std::ops::Mul<base::matrix::Matrix<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic>>>` bound

error[E0277]: cannot multiply `base::matrix::Matrix<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic>>` to `base::matrix::Matrix<N, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<N, base::dimension::Dynamic, base::dimension::Dynamic>>`
  --> src/linalg/polar.rs:83:37
   |
83 |                     Some(u.adjoint()*DMatrix::from_diagonal(&svd.singular_values)*u)
   |                                     ^ no implementation for `base::matrix::Matrix<N, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<N, base::dimension::Dynamic, base::dimension::Dynamic>> * base::matrix::Matrix<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic>>`
   |
   = help: the trait `std::ops::Mul<base::matrix::Matrix<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic>>>` is not implemented for `base::matrix::Matrix<N, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<N, base::dimension::Dynamic, base::dimension::Dynamic>>`
   = help: consider adding a `where base::matrix::Matrix<N, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<N, base::dimension::Dynamic, base::dimension::Dynamic>>: std::ops::Mul<base::matrix::Matrix<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic, base::vec_storage::VecStorage<<N as alga::general::ComplexField>::RealField, base::dimension::Dynamic, base::dimension::Dynamic>>>` bound

error: aborting due to 2 previous errors

For more information about this error, try `rustc --explain E0277`.
error: Could not compile `nalgebra`.


To learn more, run the command again with --verbose.

[sebcrozet]

Hi, thank you for this PR!

The problem here is that you are trying to multiply a matrix with elements of type N with a matrix with elements of type N::RealField. You will need to ensure the elements of DMatrix::from_diagonal(&svd.singular_values) are complex too. For example you can do DMatrix::from_diagonal(&svd.singular_values.map(|e| N::from_real(e))). This won't be the most efficient approach but we'll see about performance later.

[LucasCampos]

This seems to have fixed the problem!

Now I will write the serialize/deserialize code, some tests + benchmarks and ask for a new review.

[LucasCampos]

Just a small update. This pull request seems relatively ready for the beta version, but some things are still missing. Below I summarize what we have, what we need, and where I think we should go in the near future.

What is already implemented

1. A working (even if not complete, fast, or general) version of the decomposition
2. Benchmarks supporting it
3. Simple tests

Missing for the beta version are

1. Testing Serde: This one I am not sure how to do automatically
2. Running the tests using QuickCheck: I cannot compile them due to issue#659

What I think we should aim for in the final version

1. Accepting more general input. That is, allowing for a fixed size matrix, rather than the dynamic versions
2. Slightly change the interface of polar to allow fine-grained control of the construction of left/right Hermitian matrices
3. Write a more comprehensive set of tests

[LucasCampos]

The first alpha version seems to be ready. We have it working properly for DMatrix, as well as some tests and benchmarks, with everything compiling. From here, I think we should focus on making it possible to be used with fixed matrix formats. There, however, I will some help with the traits of the nalgebra.

This might be a good moment for a general code review.

There are some things I would like to discuss about as well, as they are on the API/user experience level. Those copy and expand on the last comment.

Slightly change the interface of polar to allow fine-grained control of the construction of left/right Hermitian matrices

The way I envision this is to copy the interface of the SVD module. On one hand, this might create a bit more cognitive load to the user, but on the other hand, it will allow for a much finer control of the operations, and I suspect a decent speedup in the common case on which only the rotation matrix R is needed.

I also can see the point of having another interface to create the polar decomposition, which accepts an pre-computed SVD. This one I am a bit less sure how useful will be.

Write more comprehensive set of tests for square matrices

If the input matrix is square, there are more tests we can do regarding the properties of each component in the decomposition. Namely, that the rotation matrices are indeed unitary.

[sebcrozet]

Sorry @LucasCampos for the very late answer. I've been sick on and off during the past two weeks. I'm all OK now.

Thank you for the work you have done so far!

The way I envision this is to copy the interface of the SVD module.

That sounds like a good idea. If the interface is the same as for the SVD, I don't think this will be too much of a cognitive load.

[sebcrozet]

From here, I think we should focus on making it possible to be used with fixed matrix formats.

Alright, so this is a difficult part. Before you begin, I suggest:

1. Commenting-out the cfg_attr for serialization. That way you won't have to worry about it at first.
2. Comment-out all the functions except fn try_new. That way you can focus on this.

Next, there are two steps: determine the trait bounds required by the struct Polar and then the additional trait bounds required by the impl<N: ComplexField> Polar<N>.

Trait bounds for struct Polar

First, each field should contain generic matrices. Meaning it should look like this:

pub struct Polar<N: ComplexField, R: Dim, C: Dim>
{
    /// The rotation matrix
    pub r: Option<MatrixMN<N, R, C>>,
    /// The left hermitian matrix (A = PR)
    pub p_l: Option<MatrixMN<N, R, R>>,
    /// The right hermitian matrix (A = RP)
    pub p_r: Option<MarixMN<N, C, C>>
}

Then, it is necessary to add trait bounds so that the compiler knows how to represent the data storage for those matrices, depending on R and C. This is the goal of the Allocator: DefaultAllocator trait bounds:

pub struct Polar<N: ComplexField, R: Dim, C: Dim>
where Allocator: DefaultAllocator<N, R, C> + DefaultAllocator<N, R, R> + DefaultAllocator<N, C, C>
{

The rule is simple: one DefaultAllocator bounds for each distinct MatrixMN, i.e., the DefaultAllocator<N, R, C> bound tells the compiler how to represent a MatrixMN<N, R, C> matrix.

The trait bounds for impl<N: ComplexField> Polar<N>

This will be the most difficult part.

First, the impl becomes:

impl<N: ComplexField, R: Dim, C: Dim> Polar<N, R, C>`
where // The new trait bounds we will talk about will go in this `where` clause.

In addition try_new becomes:

    pub fn try_new(
        matrix: MatrixMN<N, R, C>,
        eps: N::RealField,
        max_niter: usize
    ) -> Option<Self> {

Then, the goal will be to include all the trait bounds you need for every operation done on try_new. Because you are using SVD::try_new internally, you will have to copy all the trait bounds of SVD::try_new, i.e., all those. Yes, this is a bit ugly, but Rust requires us to be very explicit about trait bounds so they have to be repeated on each callers. But if we're lucky, that's all you will need to add.

If it happens that that's not the only bounds you need to add, I suggest trying to determine what operation on your code is missing some trait bounds, and see what trait bounds this specific operation requires.

Of course feel free to ask if you end up blocked with some arcanic compilation errors!

[LucasCampos]

Hello.

I want to get back working on this project now that I have a little bit more free time. It seems that in the really long interval since my last commit something changed on how one builds and runs the tests in nalgebra. Specifically, seems like the tests with quickcheck are broken.

When running cargo as

cargo test

The tests compile and run without any issues. However, when using

cargo test --feature arbitrary

compilation fails, as seen in the error log below. Are there any workarounds for this issue, or are the quickcheck tests deprecated?

I have tried both Cargo/rustc 1.44.1 and nighty.

   Compiling nalgebra v0.21.1 (/Users/lucascampos/repos/nalgebra)
error[E0277]: the trait bound `na::Complex<f64>: approx::RelativeEq` is not satisfied
   --> tests/linalg/bidiagonal.rs:59:1
    |
59  | gen_tests!(complex, RandComplex<f64>);
    | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ the trait `approx::RelativeEq` is not implemented for `na::Complex<f64>`
    | 
   ::: /Users/lucascampos/.cargo/registry/src/github.com-1ecc6299db9ec823/approx-0.3.2/src/lib.rs:257:8
    |
257 |     A: RelativeEq<B> + ?Sized,
    |        ------------- required by this bound in `approx::Relative`
    |
    = note: required because of the requirements on the impl of `approx::RelativeEq` for `na::Matrix<na::Complex<f64>, na::Dynamic, na::Dynamic, na::VecStorage<na::Complex<f64>, na::Dynamic, na::Dynamic>>`
    = note: this error originates in a macro (in Nightly builds, run with -Z macro-backtrace for more info)

    (...) About 1000 more lines like this

error[E0277]: the trait bound `na::Complex<f64>: approx::RelativeEq` is not satisfied
   --> tests/linalg/tridiagonal.rs:54:1
    |
54  | gen_tests!(complex, RandComplex<f64>);
    | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ the trait `approx::RelativeEq` is not implemented for `na::Complex<f64>`
    | 
   ::: /Users/lucascampos/.cargo/registry/src/github.com-1ecc6299db9ec823/approx-0.3.2/src/lib.rs:257:8
    |
257 |     A: RelativeEq<B> + ?Sized,
    |        ------------- required by this bound in `approx::Relative`
    |
    = note: required because of the requirements on the impl of `approx::RelativeEq` for `na::Matrix<na::Complex<f64>, na::U4, na::U4, na::ArrayStorage<na::Complex<f64>, na::U4, na::U4>>`
    = note: this error originates in a macro (in Nightly builds, run with -Z macro-backtrace for more info)

error[E0277]: the trait bound `na::Complex<f64>: approx::RelativeEq` is not satisfied
   --> tests/linalg/tridiagonal.rs:54:1
    |
54  | gen_tests!(complex, RandComplex<f64>);
    | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ the trait `approx::RelativeEq` is not implemented for `na::Complex<f64>`
    | 
   ::: /Users/lucascampos/.cargo/registry/src/github.com-1ecc6299db9ec823/approx-0.3.2/src/lib.rs:257:8
    |
257 |     A: RelativeEq<B> + ?Sized,
    |        ------------- required by this bound in `approx::Relative`
    |
    = note: required because of the requirements on the impl of `approx::RelativeEq` for `na::Matrix<na::Complex<f64>, na::U2, na::U2, na::ArrayStorage<na::Complex<f64>, na::U2, na::U2>>`
    = note: this error originates in a macro (in Nightly builds, run with -Z macro-backtrace for more info)

error: aborting due to 87 previous errors

For more information about this error, try `rustc --explain E0277`.
error: could not compile `nalgebra`.

To learn more, run the command again with --verbose.

test_log.log

[metric-space]

@LucasCampos @sebcrozet was wondering if I can take on fixing this to get this merged if the original author (@LucasCampos ) has no time ?

[LucasCampos]

I would be thrilled if you did! Just be warned that the code is kinda old, and things might have changed in the base repo.

[Andlon]

I just took a brief look over the original PR. I'm a bit confused why there are three matrices in the decomposition. Usually the polar decomposition looks like

A = R S

where R is orthogonal and S is positive semi-definite.

Anyway, as @LucasCampos said, much has changed in nalgebra since the original PR. Perhaps starting from scratch but taking inspiration and/or implementation details from the PR might be a fruitful strategy?

[LucasCampos]

@Andlon,

I'm a bit confused why there are three matrices in the decomposition.

One can perform the left or the right polar decomposition. I.e., A = R S or A = S R'. Both of these forms are quite direct to compute from the SVD. The plan was to eventually expose three different functions, one for each handedness of the decomposition and another with both.

Perhaps starting from scratch but taking inspiration and/or implementation details from the PR might be a fruitful strategy?

Agreed.

[Andlon]

I've never seen the second variant before. All applications of polar decomposition I've ever seen uses the form A = R S, just as everyone uses QR decomposition instead of LQ decomposition. Perhaps we should just stick to one, as the other form seems very niche? (Or am I wrong here?)

[LucasCampos]

I have only seen the other form being used in analytical proofs, never in numerics. For my originally intended use, the form A = R S would suffice.

I guess it makes sense to implement a single version for now and see if a concrete usage of other uses emerge in the future. For what it is worth, scipy implements both forms, so there might be some use I am not aware of.

[Andlon]

Thanks for the scipy link, that's indeed interesting.

From my perspective, it would make sense to first just implement A = RS only for square matrices (I note that scipy also can compute a generalized version for rectangular matrices). This is by far the most used version, and I suspect it covers the vast majority of applications. Then perhaps we can rather wait and see if anyone ever needs anything else.

[LucasCampos]

Completely agree. I guess the final call depends on @metric-space.

[sebcrozet]

Reading this PR again today, I’m wondering if there is any benefit in having a dedicated structure for the polar decomposition?
The dedicated Polar type introduced in this PR only builds the polar decomposition, without providing any particular operations on it. So I’m thinking that we could either:

• Not have a dedicated data-structure and simply add methods like, svd.to_polar(), matrix.polar() that return the two matrices directly.
• Or keep the struct Polar but add some operations that are commonly performed using a polar decomposition.

[LucasCampos]

The design was based on that of src/linalg/svd.rs. For the specific choice, I will simply defer to whatever @metric-space and @sebcrozet think is best.

[metric-space]

Apologies for the radio silence

Perhaps starting from scratch but taking inspiration and/or implementation details from the PR might be a fruitful strategy

Makes sense and can and will do

Not have a dedicated data-structure and simply add methods like, svd.to_polar(), matrix.polar() that return the two matrices directly

This route seems the more sensible to me, but this might change given more time spent in the trench/ more thought.
Note: Commencing work on this

Logistics question:
@sebcrozet @Andlon I did want to make sure @LucasCampos's commits stay in the history so I'm avoiding creating a fresh branch off of dev. Was wondering if that sounds alright with you?

[metric-space]

Update on this. Starting off a fresh branch. I'm fiddling around with bound traits and trying to figure out if it makes sense to bypass having to use DMatrix related ops in the original PR, but use OMatrix related ops

[metric-space]

Update on this: still working through, just a bit in wandering mode when it comes to understanding type-level ops and allocator stuff

[metric-space]

As an explicit note to everyone involved in this thread: the cleaned up PR for this was put up and has been merged #1050

Thank you all. Hopefully we can close this PR