Opened 12 years ago

Closed 6 years ago

#4932 closed defect (fixed)

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

Reported by: was Owned by: jason
Priority: major Milestone: sage-6.6
Component: linear algebra Keywords:
Cc: Merged in:
Authors: Rob Beezer Reviewers: Jeroen Demeyer
Report Upstream: N/A Work issues:
Branch: 21b9712 (Commits, GitHub, GitLab) Commit: 21b9712d4d43bb2a58d717640734b06ccffbf502
Dependencies: Stopgaps:

Status badges

Description


Attachments (1)

trac_4932-remove-solve-left-LU.patch (2.7 KB) - added by rbeezer 10 years ago.

Download all attachments as: .zip

Change History (23)

comment:1 Changed 12 years ago by jason

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.

comment:2 Changed 12 years ago by jason

  • Owner changed from was to jason
  • Status changed from new to assigned

comment:3 Changed 12 years ago by jason

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
  1. scipy.linalg.lu_factor returns a much more compact answer
  1. 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.

comment:4 follow-up: Changed 12 years ago by jason

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

comment:5 in reply to: ↑ 4 Changed 12 years ago by jason

Replying to jason:

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

comment:6 Changed 12 years ago by mabshoff

  • Milestone changed from sage-3.2.3 to sage-3.3

Better luck in 3.3.

Cheers,

Michael

comment:7 Changed 12 years ago by mabshoff

  • Milestone changed from sage-3.4 to sage-3.4.1

3.4 is for ReST tickets only.

Cheers,

Michael

comment:8 Changed 11 years ago by dagss

  • Report Upstream set to N/A
  • Status changed from new to needs_work

I'm thinking of fixing this ticket as soon as I've got time for it.

At any rate, one piece of opinion:

  • solve_left_LU should be removed (the current NotImplementedError? has served as deprecation for over a year)
  • solve_left should, for numerical types, take an algorithm argument to switch between different numerical methods for solving, where LU is one of them. solve_left should
  • By default, solve_left should use some iterative steps to ensure an error below a given treshold (and, say, fall back to QR etc.).

In addition, numerical matrices should have a set_solve_algorithm method. It is often known in the outer, calling code what kind of numerical algorithm is needed to solve a system, while some inner code might want to perform the actual call to solve_left. I.e. what algorithms perform well for a matrix is typically a property of the given matrix.

comment:9 Changed 10 years ago by rbeezer

  • Status changed from needs_work to needs_info

It appears "straight" LU decomposition in SciPy/NumPy will accept non-square matrices as input, while the LU-factor routine still requires square input.

Do we want to try to support non-square systems with the straight LU decomposition? It would mean foregoing the built-in solve routine that extens the LU-factor decomposition, and possibly a decision would have to be made about what to cache: square (very compact) or non-square (more general).

comment:10 follow-up: Changed 10 years ago by rbeezer

According to:

http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html#numpy.linalg.solve

The generic "solve" from NumPy uses LU decomposition anyway, via an LAPACK routine.

So it does not appear to me that solve_left_LU adds anything to solve_left?

Changed 10 years ago by rbeezer

comment:11 Changed 10 years ago by rbeezer

  • Authors set to Rob Beezer
  • Status changed from needs_info to needs_review

This routine needs lots of work, has been "Not Implemented" for a long time, and with two routines for solving systems available (#7852), this should be removed. Patch does just that.

comment:12 in reply to: ↑ 10 Changed 10 years ago by jason

Replying to rbeezer:

According to:

http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html#numpy.linalg.solve

The generic "solve" from NumPy uses LU decomposition anyway, via an LAPACK routine.

So it does not appear to me that solve_left_LU adds anything to solve_left?

The point behind this patch is that the LU decomposition is cached in the matrix so that multiple solves using this matrix are much faster because you don't have to compute the LU decomposition each time.

comment:13 Changed 10 years ago by jason

(I mean, the point behind this function...)

comment:14 Changed 10 years ago by rbeezer

  • Status changed from needs_review to needs_info

OK, that makes sense, I was reading the code, not the doc string. Three ideas:

  • Probably should be solve_right.
  • How about a more descriptive name, like solve_right_precomputed?
  • LU decomposition now works for rectangular matrices (#10839). I'd think it would make the most sense to cache the "LU factor" matrix just for this purpose. That's the square matrix where L and U share space and the diagonal 1's are not present.

comment:15 Changed 8 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:16 Changed 7 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:17 Changed 7 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:18 Changed 7 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:19 Changed 6 years ago by jdemeyer

  • Branch set to u/jdemeyer/ticket/4932
  • Modified changed from 08/10/14 16:51:03 to 08/10/14 16:51:03

comment:20 Changed 6 years ago by jdemeyer

  • Commit set to 21b9712d4d43bb2a58d717640734b06ccffbf502
  • Reviewers set to Jeroen Demeyer
  • Status changed from needs_info to positive_review

New commits:

21b9712Trac #4932: remove solve_left_LU for matrices

comment:21 Changed 6 years ago by jdemeyer

  • Milestone changed from sage-6.4 to sage-6.6
  • Summary changed from fix solve_left_LU for matrix_double_dense, which was totally broken forever (?) to Remove solve_left_LU for matrix_double_dense, which was totally broken forever (?)

comment:22 Changed 6 years ago by vbraun

  • Branch changed from u/jdemeyer/ticket/4932 to 21b9712d4d43bb2a58d717640734b06ccffbf502
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.