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:  sage6.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: 
Description
Attachments (1)
Change History (23)
comment:1 Changed 12 years ago by
comment:2 Changed 12 years ago by
 Owner changed from was to jason
 Status changed from new to assigned
comment:3 Changed 12 years ago by
Okay, apparently lu_factor and lu in the scipy library return two completely different things.
Some notes from some experimentation:
 scipy.linalg.lu is usually significantly slower than scipy.linalg.lu_factor
 scipy.linalg.lu_factor returns a much more compact answer
 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 (Lidentity)+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 followup: ↓ 5 Changed 12 years ago by
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
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
 Milestone changed from sage3.2.3 to sage3.3
Better luck in 3.3.
Cheers,
Michael
comment:7 Changed 12 years ago by
 Milestone changed from sage3.4 to sage3.4.1
3.4 is for ReST tickets only.
Cheers,
Michael
comment:8 Changed 11 years ago by
 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 analgorithm
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
 Status changed from needs_work to needs_info
It appears "straight" LU decomposition in SciPy/NumPy
will accept nonsquare matrices as input, while the LUfactor routine still requires square input.
Do we want to try to support nonsquare systems with the straight LU decomposition? It would mean foregoing the builtin solve routine that extens the LUfactor decomposition, and possibly a decision would have to be made about what to cache: square (very compact) or nonsquare (more general).
comment:10 followup: ↓ 12 Changed 10 years ago by
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
comment:11 Changed 10 years ago by
 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
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 tosolve_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
(I mean, the point behind this function...)
comment:14 Changed 10 years ago by
 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
 Milestone changed from sage5.11 to sage5.12
comment:16 Changed 7 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:17 Changed 7 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:18 Changed 7 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:19 Changed 6 years ago by
 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
 Commit set to 21b9712d4d43bb2a58d717640734b06ccffbf502
 Reviewers set to Jeroen Demeyer
 Status changed from needs_info to positive_review
New commits:
21b9712  Trac #4932: remove solve_left_LU for matrices

comment:21 Changed 6 years ago by
 Milestone changed from sage6.4 to sage6.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
 Branch changed from u/jdemeyer/ticket/4932 to 21b9712d4d43bb2a58d717640734b06ccffbf502
 Resolution set to fixed
 Status changed from positive_review to closed
Yep, it looks like this might have happened in the transition to numpy and probably was my fault. Here are errors:
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.