Remove solve_left_LU for matrix_double_dense, which was totally broken forever (?)

[williamstein]

Component: linear algebra

Author: Rob Beezer

Branch/Commit: 21b9712

Reviewer: Jeroen Demeyer

Issue created by migration from https://trac.sagemath.org/ticket/4932

[jasongrout]

comment:1

Yep, it looks like this might have happened in the transition to numpy and probably was my fault. Here are errors:

sage: A = matrix(RDF, 3,3, [1,2,5,7.6,2.3,1,1,2,-1]); A

[ 1.0  2.0  5.0]
[ 7.6  2.3  1.0]
[ 1.0  2.0 -1.0]
sage: b = vector(RDF,[1,2,3])
sage: A.solve_left(b)       
(-0.113695090439, 1.39018087855, -0.333333333333)
sage: A.solve_left_LU(b)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

/home/grout/<ipython console> in <module>()

/home/grout/sage/local/lib/python2.5/site-packages/sage/matrix/matrix_double_dense.so in sage.matrix.matrix_double_dense.Matrix_double_dense.solve_left_LU (sage/matrix/matrix_double_dense.c:4930)()

TypeError: unsupported operand type(s) for *: 'NoneType' and 'NoneType'
sage: A.LU()

([0.0 1.0 0.0]
[1.0 0.0 0.0]
[0.0 0.0 1.0],
 [           1.0            0.0            0.0]
[0.131578947368            1.0            0.0]
[0.131578947368            1.0            1.0],
 [          7.6           2.3           1.0]
[          0.0 1.69736842105 4.86842105263]
[          0.0           0.0          -6.0])
sage: A.solve_left_LU(b)
too many axes: 2 (effrank=2), expected rank=1
ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (981, 0))

---------------------------------------------------------------------------
error                                     Traceback (most recent call last)

/home/grout/<ipython console> in <module>()

/home/grout/sage/local/lib/python2.5/site-packages/sage/matrix/matrix_double_dense.so in sage.matrix.matrix_double_dense.Matrix_double_dense.solve_left_LU (sage/matrix/matrix_double_dense.c:4995)()

/home/grout/sage/local/lib/python2.5/site-packages/scipy/linalg/basic.pyc in lu_solve((lu, piv), b, trans, overwrite_b)
     48         raise ValueError, "incompatible dimensions."
     49     getrs, = get_lapack_funcs(('getrs',),(lu,b1))
---> 50     x,info = getrs(lu,piv,b1,trans=trans,overwrite_b=overwrite_b)
     51     if info==0:
     52         return x

error: failed in converting 2nd argument `piv' of clapack.dgetrs to C/Fortran array

The first error comes from the code not computing LU before using the cached LU decomposition. The second error apparently comes from a mistake in calling scipy.

[jasongrout]

comment:3

Okay, apparently lu_factor and lu in the scipy library return two completely different things.

Some notes from some experimentation:

1. scipy.linalg.lu is usually significantly slower than scipy.linalg.lu_factor

2. scipy.linalg.lu_factor returns a much more compact answer

3. scipy.linalg.lu_factor returns an answer directly suitable for lu_solve, etc.

I say we cache the lu_factor results and then build the P,L,U from these results if needed.

The strictly lower triangular part of the returned matrix is the L, the upper triangular part is the U, so that the returned matrix is (L-identity)+U.

The pivot array gives the row swaps needed. For example, the following code constructs the new order of rows based on the pivot array piv:

arr=range(size)
for i in range(size):
    arr[i],arr[piv[i]] = arr[piv[i]],arr[i]

The P matrix will then have a 1 in the (i,arr[i]) position (or the (arr[i],i) position...).

THE LU MATRIX RETURNED FROM lu_factor MIGHT BE TRANSPOSED! It was in my case, for some reason.

[jasongrout]

comment:4

However, the results from lu_factor also appear to give spurious answers, while the results from lu seem more correct numerically...

[jasongrout]

comment:5

Replying to @jasongrout:

However, the results from lu_factor also appear to give spurious answers, while the results from lu seem more correct numerically...

never mind; it seems like that comment is not true; I get the same spurious answers from lu()

[sagetrac-mabshoff]

comment:6

Better luck in 3.3.

Cheers,

Michael

[sagetrac-mabshoff]

comment:7

3.4 is for ReST tickets only.

Cheers,

Michael