Matrix Multiplication to find unknowns in a loop

Hi,

I will first write the code, and then explain:

%%%%Data

T = [1 2 3 4]'; %%%%4*1 marix
L = [1 4 9 16]'; %%%%4*1 marix
a = ones(size(T)); %%%%4*1 marix
T2 = T.^2; %%%%4*1 marix
T3 = T.^3; %%%%4*1 marix
p = [T3 T2 T a]; %%%%4*4 marix

%%The equation is a*T^3 + b*T^2 + c*T + d
%%The below matrix multiplication finds the unknown a,b,c and d where Theta_M is a, b, c, d

Theta_M1 = (p'*p)\p'*L; %%%%4*1 marix, Direct multiplication result

%%%%%%%%%%%%%%%%%%%%%%First %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%The equation is c*T + d
%%The below matrix multiplication finds the unknown c and d where Theta_M is c, d

p1 = [T a] %%%%4*2 marix

theta1 = (p1'*p1)\p1'*L; %%%%2*1 marix

%%%%%%The equation is b*T^2 + c*T + d
%%The below matrix multiplication finds the unknown b,c and d where Theta_M is b, c, d

p2 = [T2 T a]; %%%%4*3 marix

theta2 = (p2'*p2)\p2'*L; %%%%3*1 marix

%%%%%%%%The equation is a*T^3 + b*T^2 + c*T + d
%%The below matrix multiplication finds the unknown b,c and d where Theta_M is b, c, d

p3 = [T3 T2 T a]; %%%%4*4 marix

theta3 = (p3'*p3)\p3'*L; %%%%%%%same as Theta_M1

%%%%%%store theta3 in a matrix of Theta_M2 form

Theta_M2 = zeros(4,3);

Theta_M2(:,1) = [theta1; 0; 0];
Theta_M2(:,2) = [theta2; 0;];
Theta_M2(:,3) = theta3;

The above data is the sample data. This Matlab code works well for small data, but I have 1000 data points. I would like to perform matrix multiplication step by step (i.e. after the comment First) and store the value as Theta_M2.

The results are as follows:

Theta_M1 = [-1.05693231944315e-13
             1.00000000000080
            -1.81188397618826e-12
            1.15818465928896e-12]

Theta_M2 = [5.00000000000000  1.00000000000002  -1.05693231944315e-13
           -5.00000000000000   -1.01252339845814e-13   1.00000000000080
            0                  1.04805053524615e-13  -1.81188397618826e-12
            0                  0                   1.15818465928896e-12]

Does anyone know how to do this and put it in a loop?

Thanks a ton in advance!

Answer 1

[John D'Errico]

Sigh. How many times must I explain this? Why it is bad in terms of numerical analysis?

You wrote this:

Theta_M1 = (p'*p)\p'*L; %%%%4*1 marix, Direct multiplication result

In fact, this is a problem. What happens is it squares the condition number of the matrix p. That is bad because if your problem is at all ill-conditioned (as polynomials often are) then you have made a poorly conditioned problem considerably worse.

The problem is that someone taught you that. If I don't teach you the problem, then you will cluelessly go on, and then teach someone else to use the poor form that you are now happily using. This propagates, one person teaching the next how to use poor numerical methods, because none of you knew any better. It even finds its way into papers and even text books, because the authors never had anyone tell them what was wrong with what they did.

The solution is TRIVIAL. Simply use the tools that MATLAB provides, as they are intended to be used! MATLAB uses high quality methods from linear algebra to solve the problem without needing to resort to your trick (often called the normal equations).

The backslash operator can solve these linear problems directly, as I show here:

Better_Theta_M1 = p\L;

Theta_M1
Theta_M1 =
     -1.05693231944315e-13
           1.0000000000008
     -1.81188397618826e-12
      1.15818465928896e-12

Better_Theta_M1
Better_Theta_M1 =
     0
     1
     0
     0

As you can see, there is a significant difference between the two estimates. What happened is your version amplifies any noise in the problem. In fact, there is noise, down in the least significant bits of those numbers. You can see the result here, even for a very basic, easy problem. That noise got amplified up as far as 1e-12 or so. Had your problem been more significant,with non tiny noise, your solution would quickly begin to produce garbage for answers.

So learn to use backslash directly, as I showed.

Next, what you DO NOT want to start doing is creating numbered variables in a loop. That is the path to terrible code. Yes, I know that is the crux of your question. LEARN TO USE ARRAYS! If the results will be of different sizes, then learn to use cell arrays, or arrays of structs.

help cell

And please read this:

https://www.mathworks.com/matlabcentral/answers/304528-tutorial-why-variables-should-not-be-named-dynamically-eval