numerical analysis question: inverses and LU decomposition

Simple question: Given LU decomposition (with partial pivoting) A=PLU,
is there an easy way to find A inverse from PLU?
-- 
%  Randy Yates                  % "She's sweet on Wagner-I think she'd die for Beethoven.
%% Fuquay-Varina, NC            %  She love the way Puccini lays down a tune, and
%%% 919-577-9882                %  Verdi's always creepin' from her room." 
%%%% <yates@ieee.org>           % "Rockaria", *A New World Record*, ELO   
http://home.earthlink.net/~yatescr

Randy Yates skrev:
> Simple question: Given LU decomposition (with partial pivoting) A=PLU,
> is there an easy way to find A inverse from PLU?

Back substituion? 

Rune

Rune Allnor wrote:

> Randy Yates skrev:

>>Simple question: Given LU decomposition (with partial pivoting) A=PLU,
>>is there an easy way to find A inverse from PLU?

> Back substituion? 

That may be, but with finite precision arithmetic you often lose
information in the process.  Keeping it in LU form preserves
such information as long as possible.

-- glen

"Rune Allnor" <allnor@tele.ntnu.no> writes:

> Randy Yates skrev:
>> Simple question: Given LU decomposition (with partial pivoting) A=PLU,
>> is there an easy way to find A inverse from PLU?
>
> Back substituion? 

Hi Rune,

Thank you for responding. 

Back-substitution is what you do with PLU when you want to solve A*x =
b. Actually, it's forward- followed by back-substitution. First you
solve L*y = P^T*b using forward-substitution, then you solve U*x = y
using back-substitution. (Recall that for permutation matrix P, 
P^{-1} = P^T.)

In all of this, the inverse A^{-1} is never computed. I'm asking if
there is a way to compute the inverse matrix A^{-1} given the PLU
factorization of A.
-- 
%  Randy Yates                  % "Watching all the days go by...    
%% Fuquay-Varina, NC            %  Who are you and who am I?"
%%% 919-577-9882                % 'Mission (A World Record)', 
%%%% <yates@ieee.org>           % *A New World Record*, ELO
http://home.earthlink.net/~yatescr

Randy Yates skrev:
> "Rune Allnor" <allnor@tele.ntnu.no> writes:
>
> > Randy Yates skrev:
> >> Simple question: Given LU decomposition (with partial pivoting) A=PLU,
> >> is there an easy way to find A inverse from PLU?
> >
> > Back substituion?
>
> Hi Rune,
>
> Thank you for responding.
>
> Back-substitution is what you do with PLU when you want to solve A*x =
> b. Actually, it's forward- followed by back-substitution. First you
> solve L*y = P^T*b using forward-substitution, then you solve U*x = y
> using back-substitution. (Recall that for permutation matrix P,
> P^{-1} = P^T.)
>
> In all of this, the inverse A^{-1} is never computed. I'm asking if
> there is a way to compute the inverse matrix A^{-1} given the PLU
> factorization of A.

Leaving out P to keep things simple, I'd start with

A = LU => A' = U'L' => AA' = LUU'L' = I

where A' is the inverse of A, all matrixes being n by n.

What remains is to find L' and U' given L and U.

UU' = I

Start with the lower row. Clearly  U_n,n*U'_n,n = 1 (write
out the matrix product to see why), so U'_n,n = 1/U_n,n.

Next, write out the product of the second to lowest row:

U_n-1,n-1*U'_n-1,n-1 + U_n-1,n*U'_n,n = 1.

The only unknown here is U'_n-1,n since U'_n,n was found
in the first step.

And that's how the iteration continues. Solve for L' in the
same manner, this stime starting at the top.

Rune

On 2006-10-10 08:56:28 -0300, Randy Yates <yates@ieee.org> said:

> "Rune Allnor" <allnor@tele.ntnu.no> writes:
> 
>> Randy Yates skrev:
>>> Simple question: Given LU decomposition (with partial pivoting) A=PLU,
>>> is there an easy way to find A inverse from PLU?
>> 
>> Back substituion?
> 
> Hi Rune,
> 
> Thank you for responding.
> Back-substitution is what you do with PLU when you want to solve A*x =
> b. Actually, it's forward- followed by back-substitution. First you
> solve L*y = P^T*b using forward-substitution, then you solve U*x = y
> using back-substitution. (Recall that for permutation matrix P, P^{-1} = P^T.)
> 
> In all of this, the inverse A^{-1} is never computed. I'm asking if
> there is a way to compute the inverse matrix A^{-1} given the PLU
> factorization of A.

A^-1 is the X you get from A X = I.

This is just a compact way of writing down n equations.

The first is for the first column of X and has a simple right hand side, e_1.
The second is for the second column of X and has a simple right hand side, e_2.
and so on for n columns and right hand sides. Note that e_j has a single one
in posoition j.

You are just using the PLU factors to solve the equations for the columns
of A^-1. There are compact ways of describing this that obscure its simplicity.

Anytime someone asks for the inverse there is always the question of why
and whether they really just need equation solving. If you just want equation
solving why go through the intermediate of the solving for the columns of
inverse? More work and more accumulated rounding errors for the "simplicity"
of using an inverse.

What you seem to have missed is that the inverse is just a precomputed set
of n standard equations. The set is well chosen and very convenient for
writing down algebra.

Gordon Sande <g.sande@worldnet.att.net> writes:

> On 2006-10-10 08:56:28 -0300, Randy Yates <yates@ieee.org> said:
>
>> "Rune Allnor" <allnor@tele.ntnu.no> writes:
>>
>>> Randy Yates skrev:
>>>> Simple question: Given LU decomposition (with partial pivoting) A=PLU,
>>>> is there an easy way to find A inverse from PLU?
>>> Back substituion?
>> Hi Rune,
>> Thank you for responding.
>> Back-substitution is what you do with PLU when you want to solve A*x =
>> b. Actually, it's forward- followed by back-substitution. First you
>> solve L*y = P^T*b using forward-substitution, then you solve U*x = y
>> using back-substitution. (Recall that for permutation matrix P, P^{-1} = P^T.)
>> In all of this, the inverse A^{-1} is never computed. I'm asking if
>> there is a way to compute the inverse matrix A^{-1} given the PLU
>> factorization of A.
>
> A^-1 is the X you get from A X = I.
>
> This is just a compact way of writing down n equations.
>
> The first is for the first column of X and has a simple right hand side, e_1.
> The second is for the second column of X and has a simple right hand side, e_2.
> and so on for n columns and right hand sides. Note that e_j has a single one
> in posoition j.
>
> You are just using the PLU factors to solve the equations for the columns
> of A^-1. There are compact ways of describing this that obscure its simplicity.
>
> Anytime someone asks for the inverse there is always the question of why
> and whether they really just need equation solving. If you just want equation
> solving why go through the intermediate of the solving for the columns of
> inverse? More work and more accumulated rounding errors for the "simplicity"
> of using an inverse.
>
> What you seem to have missed is that the inverse is just a precomputed set
> of n standard equations. The set is well chosen and very convenient for
> writing down algebra.

Hi Gordon,

Thanks for clarifying the situation. You make some very good points. 

What you propose is that we use PLU to solve PLU X = I so that X is A^{-1}.
That is certainly one way to compute the inverse using the PLU factorization
of A. 

What I was wondering is if there wasn't something more direct and more
computationally efficient. For example, if V is a lower-triangular
matrix with non-zero off-diagonal entries along only the jth column,
v_j+1,j, v_j+2,j, ..., v_n,n, then the inverse of V is 
-v_j+1,j, -v_j+2,j, ..., -v_n,n. 

I was wondering if we can leverage that with a similar structure for
U to "quickly" compute the inverses.
-- 
%  Randy Yates                  % "I met someone who looks alot like you,
%% Fuquay-Varina, NC            %             she does the things you do, 
%%% 919-577-9882                %                     but she is an IBM."
%%%% <yates@ieee.org>           %        'Yours Truly, 2095', *Time*, ELO   
http://home.earthlink.net/~yatescr

On 2006-10-10 11:15:17 -0300, Randy Yates <yates@ieee.org> said:

> Gordon Sande <g.sande@worldnet.att.net> writes:
> 
>> On 2006-10-10 08:56:28 -0300, Randy Yates <yates@ieee.org> said:
>> 
>>> "Rune Allnor" <allnor@tele.ntnu.no> writes:
>>> 
>>>> Randy Yates skrev:
>>>>> Simple question: Given LU decomposition (with partial pivoting) A=PLU,
>>>>> is there an easy way to find A inverse from PLU?
>>>> Back substituion?
>>> Hi Rune,
>>> Thank you for responding.
>>> Back-substitution is what you do with PLU when you want to solve A*x =
>>> b. Actually, it's forward- followed by back-substitution. First you
>>> solve L*y = P^T*b using forward-substitution, then you solve U*x = y
>>> using back-substitution. (Recall that for permutation matrix P, P^{-1} = P^T.)
>>> In all of this, the inverse A^{-1} is never computed. I'm asking if
>>> there is a way to compute the inverse matrix A^{-1} given the PLU
>>> factorization of A.
>> 
>> A^-1 is the X you get from A X = I.
>> 
>> This is just a compact way of writing down n equations.
>> 
>> The first is for the first column of X and has a simple right hand side, e_1.
>> The second is for the second column of X and has a simple right hand side, e_2.
>> and so on for n columns and right hand sides. Note that e_j has a single one
>> in posoition j.
>> 
>> You are just using the PLU factors to solve the equations for the columns
>> of A^-1. There are compact ways of describing this that obscure its simplicity.
>> 
>> Anytime someone asks for the inverse there is always the question of why
>> and whether they really just need equation solving. If you just want equation
>> solving why go through the intermediate of the solving for the columns of
>> inverse? More work and more accumulated rounding errors for the "simplicity"
>> of using an inverse.
>> 
>> What you seem to have missed is that the inverse is just a precomputed set
>> of n standard equations. The set is well chosen and very convenient for
>> writing down algebra.
> 
> Hi Gordon,
> 
> Thanks for clarifying the situation. You make some very good points.
> What you propose is that we use PLU to solve PLU X = I so that X is A^{-1}.
> That is certainly one way to compute the inverse using the PLU factorization
> of A.
> What I was wondering is if there wasn't something more direct and more
> computationally efficient. For example, if V is a lower-triangular
> matrix with non-zero off-diagonal entries along only the jth column,
> v_j+1,j, v_j+2,j, ..., v_n,n, then the inverse of V is -v_j+1,j, 
> -v_j+2,j, ..., -v_n,n.
> I was wondering if we can leverage that with a similar structure for
> U to "quickly" compute the inverses.

Perhaps you missed

>> There are compact ways of describing this that obscure its simplicity.

There is little point in calculating zero values. But the description of how
to avoid them makes things obscure.

There is the famous quote that "Premature optimization is the root of all evil"
attributed to Knuth but he was just quoting earlier sources. You are jumping
to the optimization and missing the simplicity. (I believe the earlier source
was Tony Hoare but I expect others will provide the correct citation.)

Exhibiting the explicit inverses of the triangular factors is possible and
they can then be multiplied. The result will only be of use for illustration.

"Rune Allnor" <allnor@tele.ntnu.no> writes:

> Randy Yates skrev:
>> "Rune Allnor" <allnor@tele.ntnu.no> writes:
>>
>> > Randy Yates skrev:
>> >> Simple question: Given LU decomposition (with partial pivoting) A=PLU,
>> >> is there an easy way to find A inverse from PLU?
>> >
>> > Back substituion?
>>
>> Hi Rune,
>>
>> Thank you for responding.
>>
>> Back-substitution is what you do with PLU when you want to solve A*x =
>> b. Actually, it's forward- followed by back-substitution. First you
>> solve L*y = P^T*b using forward-substitution, then you solve U*x = y
>> using back-substitution. (Recall that for permutation matrix P,
>> P^{-1} = P^T.)
>>
>> In all of this, the inverse A^{-1} is never computed. I'm asking if
>> there is a way to compute the inverse matrix A^{-1} given the PLU
>> factorization of A.
>
> Leaving out P to keep things simple, I'd start with
>
> A = LU => A' = U'L' => AA' = LUU'L' = I
>
> where A' is the inverse of A, all matrixes being n by n.
>
> What remains is to find L' and U' given L and U.
>
> UU' = I
>
> Start with the lower row. Clearly  U_n,n*U'_n,n = 1 (write
> out the matrix product to see why), so U'_n,n = 1/U_n,n.
>
> Next, write out the product of the second to lowest row:
>
> U_n-1,n-1*U'_n-1,n-1 + U_n-1,n*U'_n,n = 1.
>
> The only unknown here is U'_n-1,n since U'_n,n was found
> in the first step.
>
> And that's how the iteration continues. Solve for L' in the
> same manner, this stime starting at the top.
>
> Rune

Yes, that's exactly the sort of thing I was looking for, Rune. Thanks.
-- 
%  Randy Yates                  % "She's sweet on Wagner-I think she'd die for Beethoven.
%% Fuquay-Varina, NC            %  She love the way Puccini lays down a tune, and
%%% 919-577-9882                %  Verdi's always creepin' from her room." 
%%%% <yates@ieee.org>           % "Rockaria", *A New World Record*, ELO   
http://home.earthlink.net/~yatescr

Randy Yates skrev:

> What I was wondering is if there wasn't something more direct and more
> computationally efficient. For example, if V is a lower-triangular
> matrix with non-zero off-diagonal entries along only the jth column,
> v_j+1,j, v_j+2,j, ..., v_n,n, then the inverse of V is
> -v_j+1,j, -v_j+2,j, ..., -v_n,n.
>
> I was wondering if we can leverage that with a similar structure for
> U to "quickly" compute the inverses.

Back substitution, as I sketched it in a different post, does that with

O(n) complexity. The hard part in matrix inversion is to compute
L and U, both what numerical precision and computational complexity
is concerned. LU factorization is, I believe, of O(n^3) complexity.

Rune