Issue 222 in mpmath: QR / LQ algorithm doesn't exist in mpmath

[mpm...@googlecode.com]

Status: New
Owner: ----

New issue 222 by ken.al...@sbcglobal.net: QR / LQ algorithm doesn't exist
in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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.

Attachments:
QR.py 1.8 KB

[mpm...@googlecode.com]

Comment #1 on issue 222 by fredrik....@gmail.com: QR / LQ algorithm doesn't

exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

Thanks. Could you add some basic test code?

Does it work with complex matrices?

[mpm...@googlecode.com]

Comment #2 on issue 222 by ken.al...@sbcglobal.net: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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.

[mpm...@googlecode.com]

Comment #3 on issue 222 by ken.al...@sbcglobal.net: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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?

[mpm...@googlecode.com]

Comment #4 on issue 222 by fredrik....@gmail.com: QR / LQ algorithm doesn't

exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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.

[mpm...@googlecode.com]

Comment #5 on issue 222 by ken.al...@sbcglobal.net: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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

Attachments:
BasicTestQRmpmath.py 8.3 KB

[mpm...@googlecode.com]

Comment #6 on issue 222 by fredrik....@gmail.com: QR / LQ algorithm doesn't

exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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.

[mpm...@googlecode.com]

Comment #7 on issue 222 by ken.al...@sbcglobal.net: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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.

[mpm...@googlecode.com]

Comment #8 on issue 222 by fredrik....@gmail.com: QR / LQ algorithm doesn't

exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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.

[mpm...@googlecode.com]

Comment #9 on issue 222 by ken.al...@sbcglobal.net: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

For doctest, you mean something like this?

mp.pretty = True
A = matrix([[1, 2], [3, 4], [1, 1]])
A
Q, R = QR(mp, A, 'full', 5)
Q
R
Q * R
Q.T * Q
Q * Q.T
A = matrix([[1+0j, 2-3j], [3+j, 4+5j]])
A
Q, R = QR(mp, A)
Q
R
Q * R
Q.T * Q.conjugate()
Q.conjugate() * Q.T

[mpm...@googlecode.com]

Comment #10 on issue 222 by ken.al...@sbcglobal.net: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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

[mpm...@googlecode.com]

Comment #11 on issue 222 by ken.al...@sbcglobal.net: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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

Attachments:
SimpleTestQRmpmath.py 2.4 KB

[mpm...@googlecode.com]

Comment #12 on issue 222 by fredrik....@gmail.com: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

I have pushed the code here:

https://github.com/fredrik-johansson/mpmath/commit/83f136839dd3c61bc6b9ae051543d1f0b3658240

Please let me know if any changes should be made.

[mpm...@googlecode.com]

Comment #13 on issue 222 by ken.al...@sbcglobal.net: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

Second Example is incorrect. It should be:

>>> B = matrix([[1+0j, 2-3j], [3+j, 4+5j]])

>>> Q, R = QR(B)
>>> Q
[(-0.196116135138184 + 0.0j) (0.281767783787004 - 0.939225945956681j)]
[ (-0.98058067569092 + 0.0j) (-0.0563535567574008 + 0.187845189191336j)]
>>> R
[(-5.09901951359278 + 0.0j) (-4.31455497304005 - 4.31455497304005j)]
[ 0.0 (4.09502512437113 + 0.0j)]

>>> Q * R
[(1.0 + 0.0j) (2.0 - 3.0j)]

[(5.0 + 0.0j) (4.0 + 5.0j)]
>>> chop(Q.T * Q.conjugate())
[1.0 0.0]
[0.0 1.0]

[mpm...@googlecode.com]

Comment #14 on issue 222 by ken.al...@sbcglobal.net: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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

[mpm...@googlecode.com]

Comment #15 on issue 222 by fredrik....@gmail.com: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

My mistake, fixed now. Thanks!

[mpm...@googlecode.com]

Comment #16 on issue 222 by ken.al...@sbcglobal.net: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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?

[mpm...@googlecode.com]

Comment #17 on issue 222 by fredrik....@gmail.com: QR / LQ algorithm

doesn't exist in mpmath
http://code.google.com/p/mpmath/issues/detail?id=222

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.