QR / LQ algorithm doesn't exist in mpmath

GoogleCodeExporter commented 9 years ago

What steps will reproduce the problem?
1. QR decomposition algorithm doesn't exist in mpmath
2.
3.

What is the expected output? What do you see instead?
A -> Q, R or A -> L, Q

What version of the product are you using? On what operating system?
2.6 on windows XP

Please provide any additional information below.

I've written a QR alogrithm to factor A -> Q, R or A -> L, Q depending on the 
dimensions of A.  There are no restrictions on the content of A (other than m > 
1 and n > 1).  A can even be all zeros.

Original issue reported on code.google.com by ken72c...@yahoo.com on 29 Mar 2012 at 4:14

Attachments:

QR.py

GoogleCodeExporter commented 9 years ago

Thanks. Could you add some basic test code?

Does it work with complex matrices?

Original comment by fredrik....@gmail.com on 29 Mar 2012 at 7:53

GoogleCodeExporter commented 9 years ago

I've never dealt with complex before so really don't know, but I doubt it.� 
I'll 
look into it and get back to you.

Original comment by ken72c...@yahoo.com on 30 Mar 2012 at 5:30

GoogleCodeExporter commented 9 years ago

I now have a basic prototype QR function that works for real (mpf) and complex 
(mpc) numbers.  I need 2-3 weeks to clean it up for submission.  Are there any 
other relevant data types it should accept?

Original comment by ken72c...@yahoo.com on 19 Apr 2012 at 4:50

GoogleCodeExporter commented 9 years ago

That sounds fantastic! 

It would be nice if it worked with the interval and double precision contexts 
as well, but that's not as important.

Original comment by fredrik....@gmail.com on 19 Apr 2012 at 4:53

GoogleCodeExporter commented 9 years ago

Attached is test program for my proposed QR function that works for real (i.e. 
mpf) and complex (i.e. mpc) data types.  The QR function itself is in the top 
of the file.  I've tested the function on over 1M random matrices.  The largest 
matrix tested was a 500 x 500 (though it took almost an hour to complete on my 
laptop).  Let me know if there are any questions.

Ken

Original comment by ken72c...@yahoo.com on 18 May 2012 at 4:24

Attachments:

BasicTestQRmpmath.py

GoogleCodeExporter commented 9 years ago

Great! All that's needed is some very simple automatic test code, and doctests 
showing QR decomposition of one or two 2x2 matrices.

Do you want to prepare a patch for mpmath? If not, I can integrate your code 
this weekend.

Original comment by fredrik....@gmail.com on 18 May 2012 at 10:09

GoogleCodeExporter commented 9 years ago

Not sure I understand.  What I submitted was a test program in the form of a py 
file.  Just run it and the program will perform QR decomposition and 
verification on random real and complex matrices (up to size 26 x 26).

I'm afraid I don't know what a doctest is.  Can you explain?

Also, I'm not familiar with mpmath patches, so I guess asking you to integrate 
the code is the only choice.

Original comment by ken72c...@yahoo.com on 18 May 2012 at 5:29

GoogleCodeExporter commented 9 years ago

There should be some test code that does a small number of tests (so it doesn't 
take too long time) and raises an exception if the implementation is broken, as 
part of the test suite.

A doctest is an interactive example for the documentation that also serves as a 
test. 

But don't worry; I can do those.

Original comment by fredrik....@gmail.com on 18 May 2012 at 5:44

GoogleCodeExporter commented 9 years ago

[deleted comment]

GoogleCodeExporter commented 9 years ago

Or this?
>>> mp.pretty = True
>>> A = matrix([[1, 2], [3, 4], [1, 1]])
>>> A
[1.0  2.0]
[3.0  4.0]
[1.0  1.0]
>>> Q, R = QR(mp, A, 'full', 5)
>>> Q
[-0.301511344577764   0.861640436855329   0.408248290463863]
[-0.904534033733291  -0.123091490979333  -0.408248290463863]
[-0.301511344577764  -0.492365963917331   0.816496580927726]
>>> R
[-3.3166247903554  -4.52267016866645]
[             0.0  0.738548945875996]
[             0.0                0.0]
>>> Q * R
[1.0  2.0]
[3.0  4.0]
[1.0  1.0]
>>> Q.T * Q
[                 1.0    2.5539004282831e-22   8.93865149899085e-22]
[ 2.5539004282831e-22                    1.0  -6.23789553736906e-22]
[8.93865149899085e-22  -6.23789553736906e-22                    1.0]
>>> Q * Q.T
[                  1.0  -4.90659624428528e-22   -5.5431189219596e-22]
[-4.90659624428528e-22                    1.0  -9.19251424970607e-22]
[ -5.5431189219596e-22  -9.19251424970607e-22                    1.0]
>>> B = matrix([[1+0j, 2-3j], [3+j, 4+5j]])
>>> B
[(1.0 + 0.0j)  (2.0 - 3.0j)]
[(3.0 + 1.0j)  (4.0 + 5.0j)]
>>> Q, R = QR(mp, B)
>>> Q
[              (-0.301511344577764 + 0.0j)   (0.069579541056407 - 
0.950920394437562j)]
[(-0.904534033733291 - 0.301511344577764j)  (-0.115965901760678 + 
0.278318164225628j)]
>>> R
[(-3.3166247903554 + 0.0j)  (-5.72871554697751 - 2.41209075662211j)]
[                      0.0                (3.91964747951093 + 0.0j)]
>>> Q * R
[(1.0 + 0.0j)  (2.0 - 3.0j)]
[(3.0 + 1.0j)  (4.0 + 5.0j)]
>>> Q.T * Q.conjugate()
[                   (1.0 - 5.9115199013348e-18j)  (5.90595164813473e-18 - 
5.90780773358393e-18j)]
[(-1.81892920363318e-18 - 7.57887175081327e-19j)                    (1.0 - 
4.9969153764245e-18j)]
>>> Q.conjugate() * Q.T
[                 (1.0 - 2.72539551051788e-18j)  (1.09015820466429e-17 + 
4.54232585640879e-18j)]
[(1.86264629050057e-17 - 6.07594719779944e-19j)                   (1.0 - 
8.18303976724142e-18j)]

Original comment by ken72c...@yahoo.com on 18 May 2012 at 8:16

GoogleCodeExporter commented 9 years ago

Attached is a simple test program that runs in seconds and performs consistency 
checks.

Original comment by ken72c...@yahoo.com on 18 May 2012 at 8:51

Attachments:

SimpleTestQRmpmath.py

GoogleCodeExporter commented 9 years ago

I have pushed the code here:

https://github.com/fredrik-johansson/mpmath/commit/83f136839dd3c61bc6b9ae051543d
1f0b3658240

Please let me know if any changes should be made.

Original comment by fredrik....@gmail.com on 21 May 2012 at 1:56

GoogleCodeExporter commented 9 years ago

[deleted comment]

GoogleCodeExporter commented 9 years ago

Second Example is incorrect.  It should be:

>>> B = matrix([[1+0j, 2-3j], [3+0j, 4+5j]])
>>> B
[(1.0 + 0.0j)  (2.0 - 3.0j)]
[(3.0 + 0.0j)  (4.0 + 5.0j)]
>>> Q, R = QR(B)
>>> Q
[(-0.316227766016838 + 0.0j)    (0.134164078649987 - 0.939148550549912j)]
[(-0.948683298050514 + 0.0j)  (-0.0447213595499958 + 0.313049516849971j)]
>>> R
[(-3.16227766016838 + 0.0j)  (-4.42718872423573 - 3.79473319220206j)]
[                       0.0                (4.47213595499958 + 0.0j)]
>>> Q * R
[(1.0 + 0.0j)  (2.0 - 3.0j)]
[(3.0 + 0.0j)  (4.0 + 5.0j)]
>>> chop(Q.T * Q.conjugate())
[1.0  0.0]
[0.0  1.0]

Original comment by ken72c...@yahoo.com on 21 May 2012 at 6:51

GoogleCodeExporter commented 9 years ago

My mistake, fixed now. Thanks!

Original comment by fredrik....@gmail.com on 24 May 2012 at 9:32

GoogleCodeExporter commented 9 years ago

Not sure I understand what 'pushing the code' means.  Is the QR function 
available now if I download mpmath again and install it or is the function in 
some development version that will be released to the public at some future 
date?

Original comment by ken72c...@yahoo.com on 24 May 2012 at 6:37

GoogleCodeExporter commented 9 years ago

It's not available in a released version yet. You can go to 
http://github.com/fredrik-johansson/mpmath and "download this repository as a 
zip file" to get a copy of the latest development version.

Original comment by fredrik....@gmail.com on 24 May 2012 at 7:11

upiterbarg / mpmath

QR / LQ algorithm doesn't exist in mpmath #222