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)

Was it helpful?

Solution 2

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

OTHER TIPS

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 nth solution will be accessible via

X{n}    % == [g1(n) g2(n) g3(n)]
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top