On Nov 28, 2007 4:06 PM, Nathan Bell <wnbell@gmail.com> wrote:
On Nov 27, 2007 4:35 PM, Anand Patil <anand.prabhakar.patil@gmail.com> wrote:
That and the analogous forward multiplication with a triangular matrix are the two things I really need right now, yeah. The backend would be simple, yet painful; upper/lower, transpose and side arguments, in addition to accounting for all the possible combinations of matrix formats, and the strides when 'B' is dense, would make for a huge number of bug opportunities. That's what made me think wrapping a library would be less work than writing even 'two' routines from scratch.
I guess I don't fully understand what you need here. What I imagined was two methods for forward/backward solves with lower/upper triangular matrices in CSR format.
I was thinking the triangular matrices could be any format, and the 'b' matrix/vector could be any format, in which case the multiplicity of functions rivals that of the ordinary multiplications involving sparse matrices. BTW, for such matrices, isn't this is equivalent to Gauss-Seidel[1]
sweep in the appropriate direction. If so, then we might just move something like [2] to scipy.linsolve.
I can't see any difference... it would be nice to have the L3 version too, though, and in that case b could be in any format. Also, are you sure the sparse LU solvers wouldn't do this for you anyway?
Nope, will have to look more closely at those.
About a year ago I implemented sparsetools from scratch. I looked at the same packages you listed before and came to the conclusion that the amount of code needed was rather small and that by writing the library to fit well with SciPy was worth the effort.
Ultimately, I think making the backend code amenable to SWIG and scipy.sparse saved more effort than starting with an existing library.
OK, if you've done the research already and come to that conclusion I should just do it! Regardless of the backend (LU, gauss-seidel or from scratch) here's what I'm proposing for the Python interface: B=trimult(A, B, uplo='U', side='L', inplace=True) B=trisolve(A, B, uplo, side, inplace) A is triangular, B is rectangular or a vector. We could trim the uplo and side arguments if there were 'triangular matrix' classes that kept track of whether they're upper or lower, but that would be more trouble than it's worth, no? I'm thinking if inplace=True, B will only be overwritten if it's possible to recycle the memory efficiently, so the only safe usage is to use the return value. How does this sound? Anand
On Nov 28, 2007 5:55 PM, Anand Patil <anand.prabhakar.patil@gmail.com> wrote:
OK, if you've done the research already and come to that conclusion I should just do it! Regardless of the backend (LU, gauss-seidel or from scratch) here's what I'm proposing for the Python interface:
B=trimult(A, B, uplo='U', side='L', inplace=True) B=trisolve(A, B, uplo, side, inplace)
I just sat down to write the 'trimult' routines and realized there's no need, the sparse general matrix-matrix multiplication saves the work anyway! Sorry for being slow. I'll just work on the 'trisolve' one.
Hi all, Sorry it took me so long to get to this, not even sure if it's needed any longer, but I think the attached satisfies my earlier need for a sparse triangular solver. Here's my thinking: Gaussian elimination is trivial for sparse upper triangular matrices, so I just wrote a wrapper for scipy.linsolve.splu that makes sure splu gets an upper triangular matrix by transposing the input matrix if necessary. Does this make sense? Cheers, Anand --Begin-- # Author: Anand Patil, 2007 # License: SciPy-compatible import scipy.sparse as sp import scipy.linsolve as ls def reverse_trans(trans): if trans=='N': return 'T' else: return 'N' class sparse_trisolver(object): """__call__(b, trans='N' or 'T')""" def __init__(self, chol_fac, uplo='U'): self.uplo = uplo if uplo=='U': self.splu = ls.splu(chol_fac) else: self.splu = ls.splu(chol_fac.T) def __call__(self, b, trans='N'): if self.uplo=='L': trans = reverse_trans(trans) if len(b.shape)==1: out = self.splu.solve(b, trans) else: out = b.copy() for i in xrange(b.shape[1]): out[:,i] = self.splu.solve(b[:,i], trans) return out # Test if __name__=='__main__': from numpy import asmatrix from numpy.random import normal from numpy.linalg import cholesky, solve A_dense = asmatrix(normal(size=(5,5))) A_dense = cholesky(A_dense.T*A_dense).T A_csc = sp.csc_matrix(A_dense) A_csr = sp.csr_matrix(A_dense) U_csc = sparse_trisolver(A_csc) U_csr = sparse_trisolver(A_csr) L_csr = sparse_trisolver(A_csc.T, 'L') L_csc = sparse_trisolver(A_csr.T, 'L') B = normal(size=(5,30)) U_real = solve(A_dense, B) L_real = solve(A_dense.T, B) U_csc_sol = U_csc(B) U_csr_sol = U_csr(B) L_csc_sol = L_csc(B) L_csr_sol = L_csr(B) for upper in U_csc_sol, U_csr_sol: print abs(upper - U_real).max() for lower in L_csc_sol, L_csr_sol: print abs(lower - L_real).max() --End-- On Nov 28, 2007 7:01 PM, Anand Patil <anand.prabhakar.patil@gmail.com> wrote:
On Nov 28, 2007 5:55 PM, Anand Patil <anand.prabhakar.patil@gmail.com> wrote:
OK, if you've done the research already and come to that conclusion I
should just do it! Regardless of the backend (LU, gauss-seidel or from scratch) here's what I'm proposing for the Python interface:
B=trimult(A, B, uplo='U', side='L', inplace=True) B=trisolve(A, B, uplo, side, inplace)
I just sat down to write the 'trimult' routines and realized there's no need, the sparse general matrix-matrix multiplication saves the work anyway! Sorry for being slow. I'll just work on the 'trisolve' one.
On Jan 15, 2008 12:51 PM, Anand Patil <anand.prabhakar.patil@gmail.com> wrote:
Here's my thinking: Gaussian elimination is trivial for sparse upper triangular matrices, so I just wrote a wrapper for scipy.linsolve.splu that makes sure splu gets an upper triangular matrix by transposing the input matrix if necessary. Does this make sense?
When you had originally mentioned sparse triangular solves, I thought you wanted something more lightweight. Are you sure that spsolve doesn't already do what you want? I would imagine that a reasonably smart LU factorization method would do little work when given an upper *or* lower triangular system. -- Nathan Bell wnbell@gmail.com http://graphics.cs.uiuc.edu/~wnbell/
When you had originally mentioned sparse triangular solves, I thought you wanted something more lightweight.
Ideally, I did... but this solution is trivial to write & maintain, and I figured it would waste no floating-point ops and only relatively few memory and logic ops.
Are you sure that spsolve doesn't already do what you want? I would imagine that a reasonably smart LU factorization method would do little work when given an upper *or* lower triangular system.
I emailed the maintainer of SuperLU just now, and she said she doubted it would do anything when asked to factorize an upper triangular system, but didn't comment on a lower triangular system. Is there any way to get the lu factors out of a factorized_lu object and check? Anand
participants (2)
-
Anand Patil -
Nathan Bell