Thread Subject: How to store factorization info when backslash, \ is used to solve matrix equation

When we want to solve Ax=b, we simply use x=A\b.
Here, \ means a lots of things. It internally tests several
algorithms to solve the given problem, chooses the best
one, and yields the solution vector x.

Now, I would like to solve the following problem: Ax=b(t),
where vector b is time-variant. Again, I can solve this
problem by using x=A\b(t) 1000 times.

However, in the case, matrix A is still time-invariant. If
I simply use x=A\b(t), Matlab will test the same matrix A
1000 times and factorize the same matrix 1000 times. I want
to avoid this situation.

When / is used first time, I want to store factorization
information abut matrix A. Then, I want to use the
information repeatedly to solve Ax=b(t).

Is there any way to collect matrix factorization
information about matrix A when we try x=A\b(t)?
Is there any better way to approach this problem using
Matlab?

Thank you for all suggestions!

Evan

"evan um" <evanum@gmail.com> wrote in message <g5j2ic$ld4
$1@fred.mathworks.com>...
> When we want to solve Ax=b, we simply use x=A\b.
> Here, \ means a lots of things. It internally tests several
> algorithms to solve the given problem, chooses the best
> one, and yields the solution vector x.
>
> Now, I would like to solve the following problem: Ax=b(t),
> where vector b is time-variant. Again, I can solve this
> problem by using x=A\b(t) 1000 times.
>
> However, in the case, matrix A is still time-invariant. If
> I simply use x=A\b(t), Matlab will test the same matrix A
> 1000 times and factorize the same matrix 1000 times. I want
> to avoid this situation.
>
> When / is used first time, I want to store factorization
> information abut matrix A. Then, I want to use the
> information repeatedly to solve Ax=b(t).
>
> Is there any way to collect matrix factorization
> information about matrix A when we try x=A\b(t)?
> Is there any better way to approach this problem using
> Matlab?
>
> Thank you for all suggestions!
>
> Evan

  What sort of matrix is your A? Is it square and non-singular? If so, you
might consider just finding its inverse once and for all and applying that fixed
inverse to all your different values, b(t). This would avoid all the repetitions
in doing A\b(t) each time.

  If A is m by n with m > n, then A\b(t) involves finding a least squares
solution. As I recall, the least squares problem can be expressed as the
solution to

 (A'*A)*x = A'*b

(I hope I have that right.) Here too, A'*A is fixed and you could compute its
inverse just once, thereafter applying this same inverse to A'*b(t) as t varies
and thereby saving some computation time.

  (I give you this advice knowing that I will undoubtedly receive a lot of flak
for so subversive a suggestion as using the inverse.)

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in
message <g5javs$n7j$1@fred.mathworks.com>...
> ..........
> If A is m by n with m > n, then A\b(t) involves finding a least squares
> solution. As I recall, the least squares problem can be expressed as the
> solution to
>
> (A'*A)*x = A'*b
>
> (I hope I have that right.) Here too, A'*A is fixed and you could compute its
> inverse just once, thereafter applying this same inverse to A'*b(t) as t varies
> and thereby saving some computation time.
> ........

  One alteration! For the m > n case I should have said, calculate

 inv(A'*A)*A'

only once and apply it to b(t) repeatedly.

Roger Stafford

Thank you for all suggestions and helps!

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
wrote in message <g5javs$n7j$1@fred.mathworks.com>...
> "evan um" <evanum@gmail.com> wrote in message <g5j2ic$ld4
> $1@fred.mathworks.com>...
> > When we want to solve Ax=b, we simply use x=A\b.
> > Here, \ means a lots of things. It internally tests several
> > algorithms to solve the given problem, chooses the best
> > one, and yields the solution vector x.
> >
> > Now, I would like to solve the following problem: Ax=b(t),
> > where vector b is time-variant. Again, I can solve this
> > problem by using x=A\b(t) 1000 times.
> >
> > However, in the case, matrix A is still time-invariant. If
> > I simply use x=A\b(t), Matlab will test the same matrix A
> > 1000 times and factorize the same matrix 1000 times. I want
> > to avoid this situation.
> >
> > When / is used first time, I want to store factorization
> > information abut matrix A. Then, I want to use the
> > information repeatedly to solve Ax=b(t).
> >
> > Is there any way to collect matrix factorization
> > information about matrix A when we try x=A\b(t)?
> > Is there any better way to approach this problem using
> > Matlab?
> >
> > Thank you for all suggestions!
> >
> > Evan
>
> What sort of matrix is your A? Is it square and
non-singular? If so, you
> might consider just finding its inverse once and for all
and applying that fixed
> inverse to all your different values, b(t). This would
avoid all the repetitions
> in doing A\b(t) each time.
>
> If A is m by n with m > n, then A\b(t) involves finding
a least squares
> solution. As I recall, the least squares problem can be
expressed as the
> solution to
>
> (A'*A)*x = A'*b
>
> (I hope I have that right.) Here too, A'*A is fixed and
you could compute its
> inverse just once, thereafter applying this same inverse
to A'*b(t) as t varies
> and thereby saving some computation time.
>
> (I give you this advice knowing that I will undoubtedly
receive a lot of flak
> for so subversive a suggestion as using the inverse.)
>
> Roger Stafford
>
>

OK ... here's the flak from the "inv police" ;-)
Never ever ever multiply by inv( ). Even multiplying by
pinv( ) has a better alternative (QR factorization).
Look for the tag "do not use inv"...

See LINFACTOR on the file exchange for help on this. You
need to factorize the matrix and then re-use the factors.
It's much faster and much more accurate.

The MathWorks could really use a "factorization object" that
captures all the details from backslash, something like:

F = mldivide_algo_but_just_give_me_the_factors_please (A) ;
x = F\b ;

which would then be syntacally as simple but much faster and
much much more accurate than

F = inv (A) ; % ack!
x = F*b ;

but in the meantime, LINFACTOR is a decent stop-gap. It
handles the full and sparse case. I don't think I does the
rectangular cases (it should do so ... maybe I'll update it
sometime to use QR as well as LU and CHOL).

Oops. I mean look for the tag "QR" or "pinv"

"Tim Davis" <davis@cise.ufl.edu> wrote in message <g5kl3v$89l
$1@fred.mathworks.com>...
> OK ... here's the flak from the "inv police" ;-)
> Never ever ever multiply by inv( ). Even multiplying by
> pinv( ) has a better alternative (QR factorization).
> Look for the tag "do not use inv"...
>
> See LINFACTOR on the file exchange for help on this. You
> need to factorize the matrix and then re-use the factors.
> It's much faster and much more accurate.
>
> The MathWorks could really use a "factorization object" that
> captures all the details from backslash, something like:
>
> F = mldivide_algo_but_just_give_me_the_factors_please (A) ;
> x = F\b ;
>
> which would then be syntacally as simple but much faster and
> much much more accurate than
>
> F = inv (A) ; % ack!
> x = F*b ;
>
> but in the meantime, LINFACTOR is a decent stop-gap. It
> handles the full and sparse case. I don't think I does the
> rectangular cases (it should do so ... maybe I'll update it
> sometime to use QR as well as LU and CHOL).

Hello Tim,

  Okay, I am ready to accept your criticism, but I need ask you a specific
question first. I'll assume here that matrix A is square and is not singular.

  First of all, we do know that if we are to find a solution vector x for all
possible vectors b by matrix multiplication of the form x = C*b, then we have
absolutely no choice as to matrix C. It must be the matrix inverse of A.
Furthermore I have to assume that matlab's 'inv' function produces the most
accurate possible approximation to this inverse and does so with reasonable
efficiency.

  I believe you stated that if a factorization object were properly stored, as in
your 'linfactor', that the execution of what you loosely referred to as "F\b", an
operation which I assume is carried out by linfactor(F,b), would be "much
faster and much much more accurate". Leaving completely aside the time it
requires to compute an inverse for A and for its factorization object, do a
thousand calls on linfactor(F,b) require less processing time than a thousand
matrix multiplications C*b? Let's assume for the sake of discussion that we
are dealing with reasonably well-conditioned, full, average run-of-the-mill
matrices for A and a wide variety of possible vectors b.

  If your answer remains that the use of stored factorizations is
unquestionably faster than simple matrix multiplication, stemming from a
wide variety of A matrices and b vectors, then I stand corrected. Remember,
we are not counting the time it takes to compute either the inverse or the
factorizations themselves.

Roger Stafford

"evan um" <evanum@gmail.com> wrote in message
news:g5j2ic$ld4$1@fred.mathworks.com...
> When we want to solve Ax=b, we simply use x=A\b.
> Here, \ means a lots of things. It internally tests several
> algorithms to solve the given problem, chooses the best
> one, and yields the solution vector x.
>
> Now, I would like to solve the following problem: Ax=b(t),
> where vector b is time-variant. Again, I can solve this
> problem by using x=A\b(t) 1000 times.
>
> However, in the case, matrix A is still time-invariant. If
> I simply use x=A\b(t), Matlab will test the same matrix A
> 1000 times and factorize the same matrix 1000 times. I want
> to avoid this situation.
>
> When / is used first time, I want to store factorization
> information abut matrix A. Then, I want to use the
> information repeatedly to solve Ax=b(t).
>
> Is there any way to collect matrix factorization
> information about matrix A when we try x=A\b(t)?
> Is there any better way to approach this problem using
> Matlab?

Do you need the result x at each time step immediately (i.e. to compute
b(t+delta_t), the b for the next time step)? If you don't, create a matrix
B where each column of B is the b vector at some time. Then solve the
system X = A\B and extract the columns.

Does the A matrix have any of the properties that the LINSOLVE function
recognizes? It wouldn't allow you to skip performing the calculations
needed by the appropriate algorithm inside backslash, but it could save you
some time testing the other methods first.

How large is your A matrix? For instance, if your A matrix is small (say
100 by 100) then each factorization is likely quick and the time you may
save by not factoring the matrix multiple times may be negligible.

Is solving the system the bottleneck in your code? If solving the system
1000 times takes a minute, but some other chunk of the code takes 5 minutes,
I'd probably start by optimizing the other chunk of code.

--
Steve Lord
slord@mathworks.com

"Tim Davis" <davis@cise.ufl.edu> wrote in message
<g5klmm$d3l$1@fred.mathworks.com>...
> Oops. I mean look for the tag "QR" or "pinv"

Dear Prof. Davis,

Thank you very much for your suggestions.
I am reading your backslash book this summer!

Best,
Evan

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in
message <g5liub$luj$1@fred.mathworks.com>...
> Hello Tim,
> Okay, I am ready to accept your criticism, but I need ask you a specific
> .........

Hello again Tim,

  Before you answer the question I posed a little while ago, you should realize
that I consider it to be rather deceptive. I might just as well have posed the
following question instead. You have been given an arbitrary square n x n
matrix C. There is a set of 1000 vectors b of size n x 1 which you are to
multiply C by, but you don't know them in advance. Can you think of any
kind of advance operation with the elements of C to obtain some object F
which would allow faster calculation on the 1000 b's than with the 1000
separate C*b matrix multiplications, once the b's have been supplied? No fair
looking at the b's before preparing the object. It should only depend on C.

  You will note that the word 'inverse' is never mentioned above, simply
because that is not germane to the discussion. I claim the question has
nothing to do with solving linear equations. It is strictly a question of
whether matrix multiplication by a constant matrix with many different
vectors can be done more efficiently by other means than this repeated
matrix multiplication itself.

  Or perhaps it can be posed another way. Given an n x n matrix C, is there
any way of speeding up the multiplication C*B where B is an arbitrary n x
1000 matrix in such a way that any advance preparations will be independent
of B? If so, I think that could be a very valuable adjunct to matlab.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
wrote in ...

> Hello Tim,
>
> Okay, I am ready to accept your criticism, but I need
ask you a specific
> question first. I'll assume here that matrix A is square
and is not singular.
>
> First of all, we do know that if we are to find a
solution vector x for all
> possible vectors b by matrix multiplication of the form x
= C*b, then we have
> absolutely no choice as to matrix C. It must be the
matrix inverse of A.
> Furthermore I have to assume that matlab's 'inv' function
produces the most
> accurate possible approximation to this inverse and does
so with reasonable
> efficiency.
 
> I believe you stated that if a factorization object were
properly stored, as in
> your 'linfactor', that the execution of what you loosely
referred to as "F\b", an
> operation which I assume is carried out by linfactor(F,b),
would be "much
> faster and much much more accurate". Leaving completely
aside the time it
> requires to compute an inverse for A and for its
factorization object, do a
> thousand calls on linfactor(F,b) require less processing
time than a thousand
> matrix multiplications C*b? Let's assume for the sake of
discussion that we
> are dealing with reasonably well-conditioned, full,
average run-of-the-mill
> matrices for A and a wide variety of possible vectors b.
>
> If your answer remains that the use of stored
factorizations is
> unquestionably faster than simple matrix multiplication,
stemming from a
> wide variety of A matrices and b vectors, then I stand
corrected. Remember,
> we are not counting the time it takes to compute either
the inverse or the
> factorizations themselves.
>
> Roger Stafford
>



The problem is that inv(A)*b is fundamentally less accurate
than A\b, no matter how well you can compute inv(A). The
problem of computing inv(A) is a less well-conditioned
problem. There's lots of theoretical analysis to back that
statement up.

LINFACTOR has its drawbacks (I state those in the code),
partly because it's an m-file and inv is built-in. Also,
the forward/backsolves are lower than they could be if F
were stored in an opaque object, coming from whatever
factorization method constructed it (not a MATLAB full or
sparse matrix).

Even with using MATLAB's L and U, and using linsolve (a
MATLAB function), S*b where S=inv(A) is precomputed is a
teeny bit faster than x=U\L\(P*b) but using linsolve in
place of backslash (telling MATLAB that L is lower
triangular, and not forcing it to look at L each time to
figure that out).

I did a fairly detailed study of this in a comment on
Loren's May 16, 2007, blog ... someone posted a cite to a
paper with some faulty experimental results. I point out
the details of the flaws in their paper in the comments on
that blog.

Case in point ... for an n-by-n matrix, n = 1000, that's
symmetric positive definite, and dense, the time to compute
inv(A) is higher than chol, but S*b is a teeny bit faster.
But ... it takes 100,000 (that's one hundred thousand)
right-hand sides, one at a time, for the use of the inverse
to be faster than forward/backsolve. Hardly worth it,
especially considering that using inv is less accurate.

So yes, I realize that one factors (or inverts, gulp!) the
matrix once, then uses the factors (or inverse, gulp!) many
times with the same A but different right-hand-sides. But
my quick look shows that even just looking at the time,
alone, it's just not worth it. The factorization is lots
faster, itself, than the inverse, and it takes lots of
right-hand-sides to make that up.

So ... never multiply by inv. Ever.