implement qr()

[GoogleCodeExporter]

>>> from mpmath import *
>>> a=matrix([[1,2],[3,4]])
>>> from mpmath.linalg import LU_decomp
>>> LU_decomp(a) # this will overwrite a
(matrix(
[[mpf('3.0'), mpf('4.0')],
 [mpf('0.33333333333333331'), mpf('0.66666666666666674')]]), [1])

It's there, but a high-level interface should be created.
Should lu() return L and U separately (and p) or using one matrix (like
LU_decomp)?

>>> from mpmath.linalg import householder
>>> a=matrix([[1,2],[3,4]])
>>> householder(a)
(matrix(
[[mpf('4.16227766016838'), mpf('-4.4271887242357328')],
 [mpf('3.0'), mpf('-0.63245553203367599')]]),
[mpf('-3.1622776601683795')],
[mpf('1.4000000000000006')],
[mpf('-0.63245553203367599')])
>>> print householder.__doc__

    (A|b) -> H, p, x, res

So qr() is not so trivial to implement. (The algorithm currently reads a
2x2 matrix as a 2x1 linear overdetermined system including right side.)

Original issue reported on code.google.com by Vinzent.Steinberg@gmail.com on 16 Sep 2008 at 4:12

[GoogleCodeExporter]

> Should lu() return L and U separately (and p) or using one matrix

Separately, I think.

Original comment by fredrik....@gmail.com on 4 Oct 2008 at 7:12

[GoogleCodeExporter]

lu() is implemented in r599.

Original comment by Vinzent.Steinberg@gmail.com on 5 Oct 2008 at 7:19

[GoogleCodeExporter]

For reference:

>>> A = matrix([[1, 2], [3, 4]])
>>> P, L, U = lu(A)
>>> P
matrix(
[[mpf('0.0'), mpf('1.0')],
 [mpf('1.0'), mpf('0.0')]])
>>> L
matrix(
[[mpf('1.0'), mpf('0.0')],
 [mpf('0.33333333333333331'), mpf('1.0')]])
>>> U
matrix(
[[mpf('3.0'), mpf('4.0')],
 [mpf('0.0'), mpf('0.66666666666666674')]])
>>> P*A
matrix(
[[mpf('3.0'), mpf('4.0')],
 [mpf('1.0'), mpf('2.0')]])
>>> L*U
matrix(
[[mpf('3.0'), mpf('4.0')],
 [mpf('1.0'), mpf('2.0')]])

I won't implement qr() anytime soon, as it's not very useful and somewhat 
complicated
to do.

Original comment by Vinzent.Steinberg@gmail.com on 6 Oct 2008 at 3:47

• Changed title: implement qr()
• Added labels: Priority-Low
• Removed labels: Priority-Medium

[GoogleCodeExporter]

What about the QR algorithm for eigenvalues? (Maybe it doesn't require the 
explicit
QR factorization; I haven't looked at the algorithm in detail.)

Original comment by fredrik....@gmail.com on 6 Oct 2008 at 6:32

[GoogleCodeExporter]

I don't know much about eigenvalues, so please don't ask me. :)

Original comment by Vinzent.Steinberg@gmail.com on 7 Oct 2008 at 3:40

[GoogleCodeExporter]

Spent some work on implementing qr(). It nearly works, but I still have to hunt 
down
some bugs...

Original comment by Vinzent.Steinberg@gmail.com on 9 Apr 2009 at 1:46

• Changed state: Started

[GoogleCodeExporter]

Cool. On a related note, does the QR decomposition support complex numbers 
correctly?
(I think it needs to work with the conjugate transpose and possibly some other
modifications visavis the real version.) This should be needed for complex 
eigenvalue
iteration and least squares solving with complex coefficients.

Original comment by fredrik....@gmail.com on 9 Apr 2009 at 1:54

[GoogleCodeExporter]

Least squares solving does not seem to work with complex numbers:

>>> A = matrix([[1j,2,4],[4j-2,3,6],[1-3j,5,7j],[1,2,3]])                      
>>> b = matrix([1+1j, 10-3j, -2j, 0])
>>> lu_solve(A, b)
------------------------------------------------------------
Traceback (most recent call last):
  File "<ipython console>", line 1, in <module>
  File "mpmath/settings.py", line 135, in g
    return f(*args, **kwargs)
  File "mpmath/linalg.py", line 210, in lu_solve
    return cholesky_solve(A, b)
  File "mpmath/settings.py", line 135, in g
    return f(*args, **kwargs)
  File "mpmath/linalg.py", line 462, in cholesky_solve
    L = cholesky(A)
  File "mpmath/linalg.py", line 436, in cholesky
    if s < eps:
  File "mpmath/mptypes.py", line 225, in __gt__
    def __gt__(s, t): return s._cmp(t, mpf_gt)
  File "mpmath/mptypes.py", line 221, in _cmp
    return func(s._mpf_, t)
  File "mpmath/libmpf.py", line 523, in mpf_gt
    return mpf_cmp(s, t) > 0
  File "mpmath/libmpf.py", line 472, in mpf_cmp
    tsign, tman, texp, tbc = t
TypeError: 'mpc' object is not iterable
>>> qr_solve(A, b)
(matrix(
[['(-4.93165217326787 + 1.89082045265996j)'],
 ['(3.94442385746051 + 0.749421140800992j)'],
 ['(-2.33395815616434 + 2.23431278511554j)']]), mpf('5.2215902676785161'))
>>> x = _[0]
>>> norm_p(residual(A,x,b),2)
mpf('18.404477839694874')

cholesky() needs to be fixed for complex numbers, if possible.

Not sure whether the least square method (used by lu_solve) works for complex 
numbers
at all. Maybe A.T*A*x == A.T*b has to be changed to A.H*A*x == A.H*b.

householder() does not calculate Q at all. You can use it for QR decomposition, 
but
currently it is quite inefficient for this, because it does not ignore b 
(right-hand
side) and always solves the system, which is not needed for decomposition only. 
This
can be easily fixed using a flag.

householder() uses sign() somewhere to avoid cancellation errors, I did not 
check
whether this works for complex numbers.

Original comment by Vinzent.Steinberg@gmail.com on 9 Apr 2009 at 2:42