Tridiagonal matrix - repeated solution
in [Matlab]
|
Prev: Data Brush Tool
Next: how to add sine to sound signal
From: Roland Kruse on 31 Mar 2010 09:21 Hello, my problem involves solving a system of equations Ax=b for different RHS vectors b. A is tridiagonal / complex / sparse. I can obviously repeatedly use x = A\b, however, this does not seem to be very fast. I know that Matlab can solve this problem very efficiently if I put all vectors b in a matrix, let's call it B, and use X = A\B. Unfortunately, the RHS vectors are calculated in an iterative fashion (from x) and, consequently, not all available at the same time. I tried LU decomposition by [L,U,P,Q,R] = lu(A), x = Q*(U\(L\(P*(R\b)))) and that was somewhat faster (maybe two to three times), but I wonder if there's a faster solution for this special case ? Thanks
From: Bruno Luong on 31 Mar 2010 09:42 "Roland Kruse" <roland.kruse(a)uni-oldenburg.de> wrote in message <hovi7u$7b7$1(a)fred.mathworks.com>... > Hello, > > my problem involves solving a system of equations Ax=b for different RHS vectors b. > A is tridiagonal / complex / sparse. > I can obviously repeatedly use x = A\b, however, this does not seem to be very fast. > I know that Matlab can solve this problem very efficiently if I put all vectors b in a matrix, let's call it B, and use X = A\B. > Unfortunately, the RHS vectors are calculated in an iterative fashion (from x) and, consequently, not all available at the same time. > I tried LU decomposition by > [L,U,P,Q,R] = lu(A), > x = Q*(U\(L\(P*(R\b)))) > and that was somewhat faster (maybe two to three times), but I wonder if there's a faster solution for this special case ? > > Thanks > > Check out Tim Davis's tool: http://www.mathworks.com/matlabcentral/fileexchange/24119-dont-let-that-inv-go-past-your-eyes-to-solve-that-system-factorize Bruno
From: John D'Errico on 31 Mar 2010 09:55 "Roland Kruse" <roland.kruse(a)uni-oldenburg.de> wrote in message <hovi7u$7b7$1(a)fred.mathworks.com>... > Hello, > > my problem involves solving a system of equations Ax=b for different RHS vectors b. > A is tridiagonal / complex / sparse. > I can obviously repeatedly use x = A\b, however, this does not seem to be very fast. > I know that Matlab can solve this problem very efficiently if I put all vectors b in a matrix, let's call it B, and use X = A\B. > Unfortunately, the RHS vectors are calculated in an iterative fashion (from x) and, consequently, not all available at the same time. > I tried LU decomposition by > [L,U,P,Q,R] = lu(A), > x = Q*(U\(L\(P*(R\b)))) > and that was somewhat faster (maybe two to three times), but I wonder if there's a faster solution for this special case ? > You say that A is sparse, BUT, is A stored in sparse form? If not, then WHY NOT? If it is, then forming the factorization in advance is as good as you can do unless you are willing to write specialized, OPTIMIZED C code as a mex function. In fact, it is indeed quite fast. Of course, it matters not how fast it is, since someone will always decide it is not fast enough for their purposes. n = 10000; A = spdiags(rand(n,3),-1:1,n,n); y = rand(n,1); tic,[L,U,P,Q,R] = lu(A);toc Elapsed time is 0.020990 seconds. tic, x = Q*(U\(L\(P*(R\y))));toc Elapsed time is 0.001050 seconds. You cannot tell me that this is not fast. John
From: Roland Kruse on 31 Mar 2010 11:55 > You cannot tell me that this is not fast. > John Sure it is fast, and yes, A is SPARSE. However, I need to execute the x = Q*(U\(L\(P*(R\y)))); line for a lot of times (let's say 1E6). What actually suprised me is that simply calling x = A\b repeatedly is not that much slower than using the above approach (which, as said, is approx. 2-3 times faster for my data). In other words, Matlab solves tridiagonal systems very efficiently. Therefore I wondered if there's a specific solution for my (tridiagonal) problem: the LU decomposition is a more general approach and may consequently be not optimal.
From: Bruno Luong on 31 Mar 2010 12:07
"Roland Kruse" <roland.kruse(a)uni-oldenburg.de> wrote in message <hovr9a$97e$1(a)fred.mathworks.com>... > > You cannot tell me that this is not fast. > > John > > Sure it is fast, and yes, A is SPARSE. > However, I need to execute the > x = Q*(U\(L\(P*(R\y)))); > line for a lot of times (let's say 1E6). > > the LU decomposition is a more general approach and may consequently be not optimal. I doubt it. The most appropriated algorithm for general triangular matrix system In know is the forward/backward substitutions based on ...LU factorization. Actually the L and U are composed by 2 diagonals. You might be able to beat it by program your own mex file but I doubt you can reduce much the time compared to Matlab. Bruno |