Added Golub Kahan Bidiagonalization with Householder Reflections

[sudoWin]

References to other Issues or PRs

Brief description of what is fixed or changed

Added Golub Kahan Bidiagonalization

Other comments

#18775

Release Notes

• matrices
  • Added Golub Kahan Bidiagonalization with Householder Reflections

[sympy-bot]

✅

Hi, I am the SymPy bot (v158). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• matrices
  • Added Golub Kahan Bidiagonalization with Householder Reflections (#18797 by @sudoWin)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.6.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->


#### Brief description of what is fixed or changed
Added Golub Kahan Bidiagonalization

#### Other comments

#18775

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* matrices
  * Added Golub Kahan Bidiagonalization with Householder Reflections
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[sudoWin]

I only added to the code of matrices.py, eigen.py, test_eigen.py, and durations.json, but for some unknown reason, tests in number theory are failing.

[sudoWin]

What should I do know? Can you please help a little @sylee957 . The tests that fail have no relation with my PR.

https://travis-ci.org/sympy/sympy/builds/659763753?utm_source=github_status&utm_medium=notification

[sylee957]

I think that the test should be rerun as I deferred the issue in #18800

[sylee957]

Can the algorithm be modified as a decomposition method to return the transformation matrix U and V?

[sudoWin]

Yes, I will modify it and submit a PR soon.

[sylee957]

I have a few suggestions

1. You should remove the evalf logic inside the code because sympy is not a numeric library and converting to floating points in priori can always be done from the user-side API.
   For now, even if things getting hanging for symbolic computations, I'd rather view that as an open research problem than simply disallowing it in this way because of the above reason.
   (But you should probably document that the computation can hang for symbolic matrices)
   And I also think that it's more educationally meaningful to have a method to compute stuff like 2*2 matrix bidiagonalization with radicals.

2. You should correct the docs like U * B * V.H (I'd rather believe it's U * B * V.H) and polish the code about spacing between operators. And I think that there should be more tests than naive cases. (e.g. complex integer matrix or rectangular matrix, mutable and immutable matrix)

3. I'd suggest a dual API like bidiagonalize and bidiagonal_decomposition, and make bidiagonalize to return only B while the decomposition returns all the U, B, V.
   I'd suggest this because a lot of the redundant computations for computing the U, V can be skipped for computing only the B matrix.

[sudoWin]

I have worked on your suggestions and here is the checklist,

• Remove evalf logic

• Added in the documentation, that computation can hang for big matrices

• Corrected the docs. and also latex formatted them.(If you think, i should add more, do tell me)

• Polished the code for spacing between operators

• Complex Integer Support Added

• Complex Integer Matrix Tests Added

• Added Random Tests along with naive tests

• Added Immutable Matrix Test

• Added Rectangular Matrix Tests

• Changed the API to have two functions, bidiagonalize, and bidiagonal_decomposition

[sudoWin]

For rectangular matrices, tests become too slow, so I had doubts on whether to add them or not. So should I add them.
For immutable matrices, I was unable to create a mutable matrix out of the immutable one, Is there a way to do so other than by using Matrix() ?

[sylee957]

If you have concerns about some slowdowns, I'd recommend adding them as floating point tests rather than radicals. I also anticipate much slowdowns for stuff larger than 2*2 matrix so we would only need 2*2 matrix tests for exact radicals and other tests can be made floating points.

There is as_mutable to convert the immutable matrix to its respectable type.

[sudoWin]

• Added Immutable Matrix Suport
• Added Immutable Matrix Test
• Added Rectangular Matrix Tests

[sudoWin]

Is there anything else that needs to be done?

[oscarbenjamin]

I think this can be improved.

The test code here is very confusing because the asserts are a long way away from the calculations. It would be better to put them right next to them like:

A, B, C = M.bidiagonalize_components()
assert A * B * C == M

There can sometimes be reasons for using random tests but I don't see the merit here because:

• The randomised tests are only for 2x2 matrices in which case bidiagonalisation is fairly trivial.
• Basic cases are not properly tested

We also normally try to avoid using generic simplify in testing or library code because usually a specific simplification function is better. Why is simplify needed? I guess from an example that cancel is sufficient:

In [43]: M = Matrix([[1, 2], [3, 4]])                                                                                                                         

In [44]: A, B, C = M.bidiagonal_decomposition()                                                                                                               

In [45]: A                                                                                                                                                    
Out[45]: 
⎡               2                            -6                 ⎤
⎢    1 - ───────────────         ───────────────────────────    ⎥
⎢                 9                        ⎛         9     ⎞    ⎥
⎢        1 + ───────────         (1 - √10)⋅⎜1 + ───────────⎟    ⎥
⎢                      2                   ⎜              2⎟    ⎥
⎢            (-1 + √10)                    ⎝    (-1 + √10) ⎠    ⎥
⎢                                                               ⎥
⎢            -6                              18                 ⎥
⎢───────────────────────────  - ──────────────────────────── + 1⎥
⎢          ⎛         9     ⎞             2 ⎛         9     ⎞    ⎥
⎢(1 - √10)⋅⎜1 + ───────────⎟    (1 - √10) ⋅⎜1 + ───────────⎟    ⎥
⎢          ⎜              2⎟               ⎜              2⎟    ⎥
⎣          ⎝    (-1 + √10) ⎠               ⎝    (-1 + √10) ⎠    ⎦

In [52]: A.applyfunc(cancel)                                                                                                                                  
Out[52]: 
⎡ -(-1 + √10)     -(-33 + 6⋅√10) ⎤
⎢ ────────────    ───────────────⎥
⎢  -10 + √10        -20 + 11⋅√10 ⎥
⎢                                ⎥
⎢-(-33 + 6⋅√10)     -31 + 13⋅√10 ⎥
⎢───────────────   ───────────── ⎥
⎣  -20 + 11⋅√10    -130 + 31⋅√10 ⎦

In [53]: A.applyfunc(cancel).applyfunc(radsimp)                                                                                                               
Out[53]: 
⎡ √10   3⋅√10⎤
⎢ ───   ─────⎥
⎢  10     10 ⎥
⎢            ⎥
⎢3⋅√10  -√10 ⎥
⎢─────  ─────⎥
⎣  10     10 ⎦

In [47]: cancel(cancel(A)*cancel(B)*C) == M                                                                                                                   
Out[47]: True

Maybe it's worth calling cancel on the outputs by default if this kind of simplification is usually beneficial.

I just tried this example which seems to hang:

In [1]: M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])                                                                                                         

In [2]: M.bidiagonalize()                                                                                                                                     
^C---------------------------------------------------------------------------
KeyError 

This also fails for simple symbolic examples:

In [57]: x = Symbol('x', real=True)                                                                                                                                      

In [58]: Matrix([[1, x], [x, 1]]).bidiagonalize()                                                                                                             
---------------------------------------------------------------------------
TypeError: cannot determine truth value of Relational

[sudoWin]

Thanks @oscarbenjamin for suggesting multiple changes,

1. I added a comment that described what the variables were storing, I will submit a new PR and move the asserts right next to calculations.

2. Could you please suggest which basic cases I should add? I would be more than happy to add them.

3. I added random tests because I thought it might be beneficial, I can remove it if you say so.
   What I thought while adding them was, that every time Testing is done, the code gets checked against a totally new Matrix, and I thought that would be better.

4. The reason behind using simplify is to simplify it as much as possible so that the final value should be an integer, which it MUST be.

In the above example itself, one has to apply

cancel
radsimp
cancel
product
cancel

Further, for assert to not raise an Assertion Error, the matrices have to be exactly equal,
and I wasn't sure after how many cancellations and simplifications, the result would be exactly equal, so that's why I used simplify()

5. We can use cancel on the outputs by default, but again, there's still a question that
   After how many cancel and radsimp would the output be so simplified that the user doesn't need to invoke simplify.

6. Examples of 3x3 or larger will hang for symbolic computation, and if the input is converted to a floating-point, well it can easily handle up to 25x25. I am looking at how to do this symbolically and with better algorithms,(and I have even found one and I am in the process of implementing it). Also
   As sylee957 said

You should remove the evalf logic inside the code because sympy is not a numeric library and converting to floating points in priori can always be done from the user-side API.
For now, even if things getting hanging for symbolic computations, I'd rather view that as an open research problem than simply disallowing it in this way because of the above reason.
(But you should probably document that the computation can hang for symbolic matrices)
And I also think that it's more educationally meaningful to have a method to compute stuff like 2*2 matrix bidiagonalization with radicals.
I will try speeding it up.

If you have concerns about some slowdowns, I'd recommend adding them as floating point tests rather than radicals. I also anticipate much slowdowns for stuff larger than 22 matrix so we would only need 22 matrix tests for exact radicals and other tests can be made floating points.

7. Regarding the usage of symbols, I don't have a solution, but I would try working to solve that one.

[sudoWin]

I tried doing the

In [57]: x = Symbol('x', real=True)                                                                                                                                      

In [58]: Matrix([[1, x], [x, 1]]).bidiagonalize()     

and the entire problem breaks down to:

In [51]: x = Symbol('x', real=True)                                                                                                                                                                         

In [52]: abs(1 + sqrt(1 + x**2)) >= abs(1 - sqrt(1 + x**2))                                                                                                                                                 
Out[52]: 
   ________       │   ________    │
  ╱  2            │  ╱  2         │
╲╱  x  + 1  + 1 ≥ │╲╱  x  + 1  - 1│

now this above expression is always True. Is there any function in Sympy that can solve this inequality? If there is one, I can try using that in the householder vector. Meanwhile, I will keep looking for a workaround.

[oscarbenjamin]

Is it fundamentally necessary to know which of the two expressions is bigger in order to make the algorithm work?

This is important because it will be the ultimate limitation in the algorithm for the symbolic case.

[sudoWin]

I am not sure, but after a lot of random tests and some pen and paperwork, I found out that the else block is responsible for flipping the sign of the diagonal element corresponding to the column being operated on, to make all diagonal values positive.

if v_plus.norm() <= v_minus.norm():
   v = v_plus
else:
   v = v_minus

and if we remove that, it would still work fine, the only thing is that signs of some values, get flipped. Technically it should be correct because the householder matrix just reflects about a plane/hyperplane and it seems totally valid as the definition of a bidiagonal matrix nowhere states that the diagonal elements must be positive. However, I couldn't find any official source or research paper stating this, so I am not 100% sure. What's your opinion on this, sir?

[sylee957]

I think that v_plus and v_minus can be chosen arbitrarily, but the first zero testing may not be skipped. But I think that the check between v_plus and v_minus can be relaxed as the citation says it only matters for numeric stability.

I'm also thinking that just using the direct formula

can be better than using the Golub's construction

Because it gives a worse result for symbolic construction. (Latter is Golub's construction)

[sudoWin]

I realized and agree that v_plus and v_minus can be chosen arbitrarily, however, I also think that both of them can be eliminated and we can also eliminate the check.
I have currently submitted a PR, which does the above and since the check was the main limitation, this also adds basic support for symbolic matrices.
Currently, the new PR uses Golub's Construction, however, I am currently looking whether the direct formula can be implemented with better and improved results.