Author 
Message 
Cheng Cosine science forum Guru Wannabe
Joined: 26 May 2005
Posts: 168

Posted: Tue Jul 11, 2006 11:27 pm Post subject:
? LS expression for A*x = b



A*x = b and A is MbyN where M < N.
Its leastsquares solution is written as:
xLS = A'*inv(A*A')*b
and this is easy to check A*xLS = A*A'*inv(A*A')*b = b
But how does one derive this when one has no idea how
solution looks like? BTW, I knew this can be expressed using
SVD, which is more general. But I am just curious on how
this expression, xLS = A'*inv(A*A')*b, was derived.
Thanks,
by Cheng Cosine
Jul/11/2k6 NC 

Back to top 


Helmut Jarausch science forum beginner
Joined: 08 Jul 2005
Posts: 49

Posted: Wed Jul 12, 2006 7:26 am Post subject:
Re: ? LS expression for A*x = b



Cheng Cosine wrote:
Quote:  A*x = b and A is MbyN where M < N.
Its leastsquares solution is written as:
xLS = A'*inv(A*A')*b

that's not true, see below
Quote: 
and this is easy to check A*xLS = A*A'*inv(A*A')*b = b
But how does one derive this when one has no idea how
solution looks like? BTW, I knew this can be expressed using
SVD, which is more general. But I am just curious on how
this expression, xLS = A'*inv(A*A')*b, was derived.

The leastsquares principles says: minimize
 A*x  b  where . denotes the Euclidean norm
Equivalently one can minimize the square of it, which is
(A*xb)' * (A*xb)
Now, a necessary condition is that the 1st derivative
must vanish. The first derivative in the directory of
(an arbitrary vector) v is
2*(A*xb)' * A*v
For this to vanish for all v one must have
(take transposes and divide by 2)
A'*(A*xb) = 0 or
(A'*A)*x = A'*b or
x= inv(A'*A)*A'*b

Helmut Jarausch
Lehrstuhl fuer Numerische Mathematik
RWTH  Aachen University
D 52056 Aachen, Germany 

Back to top 


Duncan Muirhead science forum addict
Joined: 08 Oct 2005
Posts: 70

Posted: Wed Jul 12, 2006 8:39 am Post subject:
Re: ? LS expression for A*x = b



On Wed, 12 Jul 2006 09:26:58 +0200, Helmut Jarausch wrote:
Quote:  Cheng Cosine wrote:
A*x = b and A is MbyN where M < N.
Its leastsquares solution is written as:
xLS = A'*inv(A*A')*b
that's not true, see below
and this is easy to check A*xLS = A*A'*inv(A*A')*b = b
But how does one derive this when one has no idea how
solution looks like? BTW, I knew this can be expressed using
SVD, which is more general. But I am just curious on how
this expression, xLS = A'*inv(A*A')*b, was derived.
The leastsquares principles says: minimize
 A*x  b  where . denotes the Euclidean norm
Equivalently one can minimize the square of it, which is
(A*xb)' * (A*xb)
Now, a necessary condition is that the 1st derivative
must vanish. The first derivative in the directory of
(an arbitrary vector) v is
2*(A*xb)' * A*v
For this to vanish for all v one must have
(take transposes and divide by 2)
A'*(A*xb) = 0 or
(A'*A)*x = A'*b or
x= inv(A'*A)*A'*b

There's also the completingthesquare argument:
assuming A'*A invertible:
(A*xb)'*(A*xb) = x'*A'*A*x  2*x'*A'*b + b'*b
= (xP*A'b))'*A'*A*(xP*A'b)) + b'*(IA*P*A)*b
where P = inv(A'*A)
Since A'*A >=0, the minimising x is P*A'b
Duncan 

Back to top 


Cheng Cosine science forum Guru Wannabe
Joined: 26 May 2005
Posts: 168

Posted: Wed Jul 12, 2006 8:56 am Post subject:
Re: ? LS expression for A*x = b



"Helmut Jarausch" <jarausch@igpm.rwthaachen.de> wrote in message
news:4hjmi2F1ruh7jU1@news.dfncis.de...
Quote:  Cheng Cosine wrote:
A*x = b and A is MbyN where M < N.
Its leastsquares solution is written as:
xLS = A'*inv(A*A')*b
that's not true, see below
and this is easy to check A*xLS = A*A'*inv(A*A')*b = b
But how does one derive this when one has no idea how
solution looks like? BTW, I knew this can be expressed using
SVD, which is more general. But I am just curious on how
this expression, xLS = A'*inv(A*A')*b, was derived.
The leastsquares principles says: minimize
 A*x  b  where . denotes the Euclidean norm
Equivalently one can minimize the square of it, which is
(A*xb)' * (A*xb)
Now, a necessary condition is that the 1st derivative
must vanish. The first derivative in the directory of
(an arbitrary vector) v is
2*(A*xb)' * A*v
For this to vanish for all v one must have
(take transposes and divide by 2)
A'*(A*xb) = 0 or
(A'*A)*x = A'*b or
x= inv(A'*A)*A'*b

Thanks, but I asked for A*x = b and A is MbyN where M < N.
When M < N, A'*A is of NbyN, which is not fullranked, and therefore
inv(A'*A) does not exist.
by Cheng Cosine
Jul/12/2k6 NC 

Back to top 


Peter Spellucci science forum Guru
Joined: 29 Apr 2005
Posts: 702

Posted: Wed Jul 12, 2006 10:14 am Post subject:
Re: ? LS expression for A*x = b



In article <8tWsg.19451$so3.8947@southeast.rr.com>,
"Cheng Cosine" <acosine@spamfree.com> writes:
Quote: 
A*x = b and A is MbyN where M < N.
Its leastsquares solution is written as:
xLS = A'*inv(A*A')*b
and this is easy to check A*xLS = A*A'*inv(A*A')*b = b
But how does one derive this when one has no idea how
solution looks like? BTW, I knew this can be expressed using
SVD, which is more general. But I am just curious on how
this expression, xLS = A'*inv(A*A')*b, was derived.
Thanks,
by Cheng Cosine
Jul/11/2k6 NC

since A*x=b is underdetermined, you seek the minimum norm least squares solution.
we assume rank(A)=M, otherwise the problem may be unsolvable using the
approach that follows now (and of course in your formula, inv(A'*A) makes no
sense)
we write the problem as
minimize x'*x/2
subject to A*x=b (! observe that we assume compatibility of the
constraints here, hence the rank assumption)
we write down the Lagrange multiplier necessary (and in this convex regular case
sufficient) optimality conditions:
x  A'*lambda = 0 lambda are the Lagrange multipliers for the equality
constraints
A*x = b
or in matrix form
[ I A' ] [x ] = [0]
[ A O ] [lambda] [b]
now we invert the block matrix using block Gauss elimination which is possible
since I is invertible
[ I A' ] = [ I O ] [ I A' ]
[ A O ] [ A I ] [ O A*A' ]
inverting the product
inv(...) = [ I A'*(inv(A*A') ] [ I O ]
[ O inv(A*A') ] [ A I ]
=
[ I  A'*inv(A*A')*A A'*inv(A*A') ]
[ inv(A*A')*A inv(A*A') ]
and multiplying this with the right hand side
[ 0 ]
[ b ]
you get your formula for xLS.
but this is not the way this should be computed. you use the QRdecomposition of
A' :
Q*A' = [ R ; O ] and let y=Q'*x = [y1;y2] (matlab notation)
with y1 in R^M.
the y2=0, R'*y1=b is the minimum norm solution and xLS=Q*y.
here you can indeed update, adding roews to A means adding columns to A',
this is even rather easy. but as said already, in your application I see
little sense in considering this.
hth
peter 

Back to top 


Google


Back to top 



The time now is Sun Nov 18, 2018 4:03 pm  All times are GMT

