Hi everybody,

-I'm trying to program the ordinary least squares method, and I must
solve this matrix system:

e = h - A* x
x is the unknown.

After the resolution of the system and when I try the verification, x
doesn't resolve it!
Why?

The program:
% ------------------------------------------
% ------------------------------------------
clc;

M = 1; % arbitrary
N = 5; % arbitrary
P = 5; % arbitrary but P >= Q
Q = 4; % arbitrary but P >= Q

e = zeros (M* (N+ P+ 1), 1); % arbitrary = 0
h = ones (M* (N+ P+ 1), 1); % arbitrary = ones
A = eye (M* (N+ P+ 1), P+ M* (Q+ 1)); % arbitrary = eye

x = ((A' * A)^(-1)) * (A' * (h- e)),
% the size of x is (P+ M* (Q+ 1), 1);
% (A' * A) to transform into a square matrix for inversion

e == (h - (A * x)),
% Normally after this, the result is 1
% Why it doesn't!?

% ------------------------------------------
% ------------------------------------------

Thanks in advance!

In article <eeff36b.-1@webx.raydaftYaTP>, "Zaki N."
<zaki_nabil@caramail.com> wrote:

> Hi everybody,
> 
> -I'm trying to program the ordinary least squares method, and I must
> solve this matrix system:
> 
> e = h - A* x
> x is the unknown.
> 
> After the resolution of the system and when I try the verification, x
> doesn't resolve it!
> Why?
> 
> The program:
> % ------------------------------------------
> clc;
> 
> M = 1; % arbitrary
> N = 5; % arbitrary
> P = 5; % arbitrary but P >= Q
> Q = 4; % arbitrary but P >= Q
> 
> e = zeros (M* (N+ P+ 1), 1); % arbitrary = 0
> h = ones (M* (N+ P+ 1), 1); % arbitrary = ones
> A = eye (M* (N+ P+ 1), P+ M* (Q+ 1)); % arbitrary = eye
> 
> x = ((A' * A)^(-1)) * (A' * (h- e)),
> % the size of x is (P+ M* (Q+ 1), 1);
> % (A' * A) to transform into a square matrix for inversion
> 
> e == (h - (A * x)),
> % Normally after this, the result is 1
> % Why it doesn't!?
> % ------------------------------------------
> 
> Thanks in advance!
-----------
Hello Zaki,

  Two things are wrong here.  To begin with, you're requiring

e == (h - (A * x))

to be exactly true, which makes no allowance for round off errors.  You
need to change this to allow an error within a certain judiciously
selected tolerance.

  However, the main problem is that, in effect, there are eleven linear
equalities here but only ten variables in x to adjust.  The idea in least
squares is to find the value of x which produces the smallest sum of
squares in the expression A*x-(h-e).  In the above circumstances it will,
in general, be impossible to make it identically zero.

(Remove "xyzzy" and ".invalid" to send me email.)
Roger Stafford

Thanks for your reply Roger.

-The problem isn't coming from the ordinary least squares method or
not.
We are in front of a problem of resolving a matrix system:

e = h - A* x

Where "x" is the unknown, and "h", "A", and "e" can be 'anything'
known.

After the resolution of the system and when I try the verification, x
doesn't resolve it!

Where is the mistake?
Please give me another method to resolve it.

-The program:

% ------------------------------------------
% ------------------------------------------
clc;

M = 1; % arbitrary
N = 5; % arbitrary
P = 5; % arbitrary but P >= Q
Q = 4; % arbitrary but P >= Q

e = zeros (M* (N+ P+ 1), 1); % arbitrary = 0
h = ones (M* (N+ P+ 1), 1); % arbitrary = ones
A = eye (M* (N+ P+ 1), P+ M* (Q+ 1)); % arbitrary = eye

x = ((A' * A)^(-1)) * (A' * (h- e)),
% the size of x is (P+ M* (Q+ 1), 1);
% (A' * A) to transform into a square matrix for inversion

e == (h - (A * x)),
% Normally after this, the result is 1
% Why it doesn't!?

% ------------------------------------------
% ------------------------------------------

"Zaki N." <zaki_nabil@caramail.com> wrote in message 
news:eeff36b.1@webx.raydaftYaTP...
> Thanks for your reply Roger.
>
> -The problem isn't coming from the ordinary least squares method or
> not.
> We are in front of a problem of resolving a matrix system:
>
> e = h - A* x
>
> Where "x" is the unknown, and "h", "A", and "e" can be 'anything'
> known.
>
> After the resolution of the system and when I try the verification, x
> doesn't resolve it!

For there to be a solution x to the system A*x = h (i.e. the e vector is 
exactly the zero vector) you would need to be able to write h as a linear 
combination of the columns of A -- the coefficients in the linear 
combination are the elements of x.  Note that no column of A contains a 
nonzero value in the last element.  Therefore, no matter what coefficients 
you choose, the linear combination of the columns of A will always have a 
zero as the last element of the resulting vector.  Since h has a 1 as its 
last element, we can't express h as a linear combination of the columns of A 
and so the system A*x = h has no exact solution.

Since there is no solution x that satisfies A*x = h exactly, we try to solve 
it in a least-squares sense -- minimize the norm of (A*x - h) [which is -e 
in your terminology.]  In this case, the x that minimizes this norm is the 
ones vector that the code below produces.  The residual (A*x - h) is:

[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; -1]

The last element of the residual is always going to be -1, from the 
explanation above.  Therefore to minimize the norm, we minimize the 
magnitudes of the other elements of the residual -- and 0 is as small as the 
magnitudes can get.

> Where is the mistake?

There is no mistake.

> Please give me another method to resolve it.

Just a suggestion:  don't use ^(-1) to invert your matrix.  If you need to 
invert the matrix, use the INV function.  However, you don't need to invert 
the matrix -- simply use the backslash operator (\) instead.  See HELP 
MLDIVIDE for more information on the backslash operator.

-- 
Steve Lord
slord@mathworks.com

>
> -The program:
>
> % ------------------------------------------
> % ------------------------------------------
> clc;
>
> M = 1; % arbitrary
> N = 5; % arbitrary
> P = 5; % arbitrary but P >= Q
> Q = 4; % arbitrary but P >= Q
>
> e = zeros (M* (N+ P+ 1), 1); % arbitrary = 0
> h = ones (M* (N+ P+ 1), 1); % arbitrary = ones
> A = eye (M* (N+ P+ 1), P+ M* (Q+ 1)); % arbitrary = eye
>
> x = ((A' * A)^(-1)) * (A' * (h- e)),
> % the size of x is (P+ M* (Q+ 1), 1);
> % (A' * A) to transform into a square matrix for inversion
>
> e == (h - (A * x)),
> % Normally after this, the result is 1
> % Why it doesn't!?
>
> % ------------------------------------------
> % ------------------------------------------ 

Thank you!

Even if I didn't understood all your explanations!