unable to solve the following equations in MATLAB

Question

Following is a part of a code I am writing in matlab, here I want to carry out the following simple mathematical operation [A][X]=[B], where [X] is unknown. In my case I have length of k around 1600000. So all I want is to get the values for g1,g2 and g3 for each element on the array. I have tried the following

k31 = omega3./(d)      
k32 = omega3_2./(d)    
A   = [2,1,5;-2,-1,-5];    
X   = [g1;g2;g3];

for ii = 1:length(k31)    
    B = [k31(ii); k32(ii)];    
    X = pinv(A).*B;    
end

display(g1,g2,g3)

I am using pseudo-inverse so basically I can get a solution for each X and I have made some edit there....and x is unknown, MATHEMATICALLY it can be done, but am not able to code it

Also how do I plot the values of g1 g2 g3 with x and y as follows scatter(x(1:end-1), y(1:end-1), 5, g1); scatter(x(1:end-1), y(1:end-1), 5, g2) and scatter(x(1:end-1), y(1:end-1), 5, g3)

What happens, or does not happen? You don't mention the results and expected results. — , Nov 08 '12 at 10:05
You have more unknowns (3) than equations (2). You'll either need to find another equation or do a least squares fit. — duffymo, Nov 08 '12 at 10:07
This code makes no sense...`A` and the initial `X` don't get used at all...? Also, you overwrite `X` on every iteration, making the outcome of the loop essentially equal to `X = pinv(R).*[k31(end); k32(end)];`...? — Rody Oldenhuis, Nov 08 '12 at 10:13
@TimN: The least squares fit is that solution that lies closest to the row-space of the design matrix. — Rody Oldenhuis, Nov 08 '12 at 10:16
You don't understand least squares fit. It'll give you the values that represent the straight line that minimizes the mean square error. — duffymo, Nov 08 '12 at 10:18
@dufymo: That is actually only one application of least squares; LSQ is a much broader concept. And...I meant column space :) — Rody Oldenhuis, Nov 08 '12 at 10:29
I get it. You seem to be struggling, and it's only a comment, so I kept it simple. — duffymo, Nov 08 '12 at 10:30
@Rody - no, the only solution with mean square error zero is the exact, analytic solution where the number of equations and unknowns are equal. Any other solution will find it impossible to pass through every point with a single straight line, so the MSE > 0. — duffymo, Nov 08 '12 at 10:39
@dufymo: I think you're messing up who you are replying to :) There's an awfully useful script to prevent this [here](http://userscripts.org/scripts/show/68252) (or on SO Apps) — Rody Oldenhuis, Nov 08 '12 at 10:41
@duffymo You wrote: "You have more unknowns (3) than equations (2)". In such a case, there are infinitely many solutions with MSE = 0. With one unknown less (lower degree of freedom), there would be exactly one solution with MSE = 0. With only one unknown, least-squares could be used to minimise MSE. — Tim, Nov 08 '12 at 10:44
Sorry, you have no idea what you're talking about. Yes, there are an infinite number of possible straight line solutions, but all of them have higher MSE than the least squares solution. — duffymo, Nov 08 '12 at 11:25

score 1 · Answer 1 · edited Nov 08 '12 at 11:24

1

I have to make a few assumptions here, so bear with me.

I suspect you want to do this:

k31 = omega3./(d)      
k32 = omega3_2./(d)    
A   = [2,1,5;-2,-1,-5];    

X = cell(length(k31),1);
for ii = 1:length(k31)            
    X{ii} = A\[k31(ii); k32(ii)];    
end

which uses the backslash operator instead of inv or pinv. Type help slash to get more information.

The backslash operator will usually be much faster and more accurate than inv or pinv. It is also a lot more flexible -- your case is underdetermined (you are 1 equation short of being able to solve explicitly), in which case the backslash operator will find you a least-squares solution.

Note how I save all results in a cell-array X. This means that the n^th solution will be accessible via

X{n}    % == [g1(n) g2(n) g3(n)]

edited Nov 08 '12 at 11:24

Eitan T

32,660
14
72
109

answered Nov 08 '12 at 10:21

Rody Oldenhuis

37,726
7
50
96

" in which case the backslash operator will find you a least-squares solution." - exactly. – duffymo Nov 08 '12 at 10:31
@Rody Oldenhuis I dont understand. Well in the end I require an array of values of g1, g2, and g3. which I am going to use to plot some data. can you help me in getting that – Rohan Chakrabarty Nov 08 '12 at 10:32
@RohanChakrabarty: Each vector in `X` will be a 3-element least squares fit to each system `Ax = [k31;k32]`. You can convert this to a plain array by issuing `X = [X{:}]`, which will give you a `3-by-numel(k32)` matrix of solutions. The `N`th column of this matrix is the least-squares solution to the `N`th equation. – Rody Oldenhuis Nov 08 '12 at 10:39
@Rody if I am doing this A=[2,1,5;-2,-1,-5]; G=A\[5;9] I get the following G = 0 0 -0.4000 so I am not able to understand why the other two elements become zero – Rohan Chakrabarty Nov 08 '12 at 10:53
when I do the following A=[2,1,5;-2,-1,-5]; G=A\[5;9] display(G) I get G =[0;0;-0.4], so I am not able to understand why first two elements become zero – Rohan Chakrabarty Nov 08 '12 at 10:56
@RohanChakrabarty: That's because `A*G` is the closest `G` can come to `B=[5;9]`. Basically, given your `A`, you're asking Matlab to solve an impossible set of equations: `2*x+1*y+5*z=5` (first row of `A`), while `-2*x-1*y-5*z=9` (second row of `A`), which should equal `-5` according to the first row of `A`. Therefore, this is impossible, and `G==[x y z]==[0 0 -0.4];` is the best you could hope for. – Rody Oldenhuis Nov 08 '12 at 11:02
@RodyOldenhuis Fixed a funny typo. Sorry for being a grammar freak :) – Eitan T Nov 08 '12 at 11:24
@RodyOldenhuis How do I plot the values of g1 g2 g3 with x and y as follows scatter(x(1:end-1), y(1:end-1), 5, g1); scatter(x(1:end-1), y(1:end-1), 5, g2) and scatter(x(1:end-1), y(1:end-1), 5, g3) – Rohan Chakrabarty Nov 08 '12 at 12:34

julietKiloRomeo · Accepted Answer · 2012-11-08T16:27:23.063

0

In my opinion you are better off creating a pseudo inverse from the singular value decomposition. Moore-Penrose pseudo inverse will work, but I think it gives some weird results sometimes. The least squares can be unstable, especially since the rank of your example matrix is not full.

Also; don't calculate the inverse for every iteration!

Here is an example, where rank-deficiency of A is taken care of:

A   = [2,1,5;-2,-1,-5];

% pseudo inverse of A
[m,n]   = size(A);
% get SVD
[U,S,V] = svd(A);
% check rank
r       = rank(S);
SR      = S(1:r,1:r);
% make complete if rank is not full
SRc     = [SR^-1 zeros(r,m-r);zeros(n-r,r) zeros(n-r,m-r)];
% create inverse
A_inv   = V*SRc*U.';

X=[];    
for i = 1:1600
    % this actually takes most of the time
    k31 = rand(1000, 1);
    k32 = rand(1000, 1);
    B   = [k31';k32'];
    % X is overwritten on every loop...
    X   = [X A_inv*B];
end

N_xy = 1000;
x    = rand(N_xy,1);
y    = rand(N_xy,1);
g1   = X(1,1:N_xy);

figure(1), clf
scatter(x,y,5,g1)

I didn't plot all 1600000 points since that can't be what you want

edited Nov 08 '12 at 16:27

answered Nov 08 '12 at 11:46

julietKiloRomeo

178
7

thnx...let me try...I did the following R=[2,1,5;-2,-1,-5]; for ii = 1:length(k31) X{ii} = pinv(R)*[k31(ii); k32(ii)]; g1(ii)=X(1); g2(ii)=X(2); g3(ii)=X(3); end Then all I wanted to do was plot for each g1 as follows scatter(x(1:end-1), y(1:end-1), 5, g1); but I am not getting it... – Rohan Chakrabarty Nov 08 '12 at 12:09
How do I plot the values of g1 g2 g3 with x and y as follows scatter(x(1:end-1), y(1:end-1), 5, g1); – Rohan Chakrabarty Nov 08 '12 at 12:30
Well, you haven't explained what x and y are, but see my edit above. I am sure you understand that you have to replace all the rand's with whatever it is you want to do? :) – julietKiloRomeo Nov 08 '12 at 16:28

unable to solve the following equations in MATLAB

2 Answers2