Elias Kyriakides <elias_ky...@hotmail.com> wrote in message <3AC15EAA...@hotmail.com>...
Dear friends,
Does anybody know how to get the actual code for the pseudoinverse in MATLAB? I am really interested to see the whole thing plus the Singular Value Decomposition SVD) and find out how it works.
I would be really grateful if somebody has it or knows how i can get it.
Thanks in advance,
Elias
The SVD approach of pinv() is unstable.
I developed an algorithm that solves the problem. It’s based on Tkhonov regularization.
I will make the derivation available as a pdf file, it’s part of my PhD thesis . But if you need it now then send me an email: dvasilakyATgmail.com
For now, I will give two examples and just the code for you to try it out on your problem.
%Example 1:
format long;
A = [1 0;0 0;0 0];
b=[1;1;1];
pA = pinv(A);
x = pA*b
dA = [1 0;0 10e-6;0 0];
pdA = pinv(dA);
x1 = pdA*b
%Here is the code. To try it out, set your matrix to D in my code. %*************start code*********************
D=A;
cst=.000000001; %regularization parametr, similar to tol
[r,c]=size(D); %get the number of row and columns of X
Y=eye(r);
Q(:,1)=D(:,1);
I=eye(r); %generate Identity matrix of dim rxr
sum=zeros([r,r]); %initialize sum with a rxc matrix with zeros
for k=1:c-1,
sum=sum + Q(:,k)*Q(:,k)'/(cst + D(:,k)'*Q(:,k));
Q(:,k+1)= (I-sum)*D(:,k+1);
end
for j=1:r,
S=zeros([r,1]); M(c,j)=((Q(:,c)'*Y(:,j)))/(cst + D(:,c)'*Q(:,c));
for i=1:c-1,
S=S + M(c-i+1,j)*D(:,c-i+1);
dm= 1/(cst + D(:,c-i)'*Q(:,c-i));
M(c-i,j)=(dm*Q(:,c-i)')*(Y(:,j) - S);
end
end
%*************end code*******
x3 = M*b
%Example 2 ill-conditioned Hilbert matrix
%OK A was rank deficient, how about ill-conditioned example where MatLab pinv() gives
%nonsense but my algorithm works.
% use pinv() to solve Dx = b and Dx = b1 where b1=b+db;
%then use my algorithm to solve Dx=b1 and look at the results.
%Please try this:
format long;
D=hilb(12);
x=ones(12,1); % create a sulution of 1s.
b=D*x; % b is the RHS and we know the solution x
db=[ %perturb b by db
.00000000000001;
.00000000000001;
.00000000000002;
-.0000000000001;
.00000000000001;
.00000000000001;
.00000000000002;
-.0000000000001;
.00000000000001;
.00000000000001;
.00000000000002;
-.0000000000001;
];
b1=b+db;
cst=.000000001; %regularization parametr, similar to tol
[r,c]=size(D); %get the number of row and columns of X
Y=eye(r);
Q(:,1)=D(:,1);
I=eye(r); %generate Identity matrix of dim rxr
sum=zeros([r,r]); %initialize sum with a rxc matrix with zeros
for k=1:c-1,
sum=sum + Q(:,k)*Q(:,k)'/(cst + D(:,k)'*Q(:,k));
Q(:,k+1)= (I-sum)*D(:,k+1);
end
for j=1:r,
S=zeros([r,1]); M(c,j)=((Q(:,c)'*Y(:,j)))/(cst + D(:,c)'*Q(:,c));
for i=1:c-1,
S=S + M(c-i+1,j)*D(:,c-i+1);
dm= 1/(cst + D(:,c-i)'*Q(:,c-i));
M(c-i,j)=(dm*Q(:,c-i)')*(Y(:,j) - S);
end
end
x4 =M*b1
Elias Kyriakides <elias_ky...@hotmail.com> wrote in message <3AC15EAA...@hotmail.com>...
Dear friends,
Does anybody know how to get the actual code for the pseudoinverse in MATLAB? I am really interested to see the whole thing plus the Singular Value Decomposition SVD) and find out how it works.
I would be really grateful if somebody has it or knows how i can get it.
Thanks in advance,
Elias
The SVD approach of pinv() is unstable.
I developed an algorithm that solves the problem. It’s based on Tkhonov regularization.
I will make the derivation available as a pdf file, it’s part of my PhD thesis . But if you need it now then send me an email: dvasilakyATgmail.com
For now, I will give two examples and just the code for you to try it out on your problem.
%Example 1:Daniel Vasilaky posted a photo at flickr.com
format long;
A = [1 0;0 0;0 0];
b=[1;1;1];
pA = pinv(A);
x = pA*b
dA = [1 0;0 10e-6;0 0];
pdA = pinv(dA);
x1 = pdA*b
%Here is the code. To try it out, set your matrix to D in my code. %*************start code*********************
D=A;
cst=.000000001; %regularization parametr, similar to tol
[r,c]=size(D); %get the number of row and columns of X
Y=eye(r);
Q(:,1)=D(:,1);
I=eye(r); %generate Identity matrix of dim rxr
sum=zeros([r,r]); %initialize sum with a rxc matrix with zeros
for k=1:c-1,
sum=sum + Q(:,k)*Q(:,k)'/(cst + D(:,k)'*Q(:,k));
Q(:,k+1)= (I-sum)*D(:,k+1);
end
for j=1:r,
S=zeros([r,1]); M(c,j)=((Q(:,c)'*Y(:,j)))/(cst + D(:,c)'*Q(:,c));
for i=1:c-1,
S=S + M(c-i+1,j)*D(:,c-i+1);
dm= 1/(cst + D(:,c-i)'*Q(:,c-i));
M(c-i,j)=(dm*Q(:,c-i)')*(Y(:,j) - S);
end
end
%*************end code*******
x3 = M*b
%Example 2 ill-conditioned Hilbert matrix
%OK A was rank deficient, how about ill-conditioned example where MatLab pinv() gives
%nonsense but my algorithm works.
% use pinv() to solve Dx = b and Dx = b1 where b1=b+db;
%then use my algorithm to solve Dx=b1 and look at the results.
%Please try this:
format long;
D=hilb(12);
x=ones(12,1); % create a sulution of 1s.
b=D*x; % b is the RHS and we know the solution x
db=[ %perturb b by db
.00000000000001;
.00000000000001;
.00000000000002;
-.0000000000001;
.00000000000001;
.00000000000001;
.00000000000002;
-.0000000000001;
.00000000000001;
.00000000000001;
.00000000000002;
-.0000000000001;
];
b1=b+db;
cst=.000000001; %regularization parametr, similar to tol
[r,c]=size(D); %get the number of row and columns of X
Y=eye(r);
Q(:,1)=D(:,1);
I=eye(r); %generate Identity matrix of dim rxr
sum=zeros([r,r]); %initialize sum with a rxc matrix with zeros
for k=1:c-1,
sum=sum + Q(:,k)*Q(:,k)'/(cst + D(:,k)'*Q(:,k));
Q(:,k+1)= (I-sum)*D(:,k+1);
end
for j=1:r,
S=zeros([r,1]); M(c,j)=((Q(:,c)'*Y(:,j)))/(cst + D(:,c)'*Q(:,c));
for i=1:c-1,
S=S + M(c-i+1,j)*D(:,c-i+1);
dm= 1/(cst + D(:,c-i)'*Q(:,c-i));
M(c-i,j)=(dm*Q(:,c-i)')*(Y(:,j) - S);
end
end
x4 =M*b1
Dear friends,Daniel Vasilaky posted a photo at flickr.com
Does anybody know how to get the actual code for the pseudoinverse in
MATLAB? I am really interested to see the whole thing plus the Singular
Value Decomposition SVD) and find out how it works.
I would be really grateful if somebody has it or knows how i can get it. Thanks in advance,
Elias
| Sysop: | Keyop |
|---|---|
| Location: | Huddersfield, West Yorkshire, UK |
| Users: | 296 |
| Nodes: | 16 (2 / 14) |
| Uptime: | 81:24:02 |
| Calls: | 6,658 |
| Calls today: | 4 |
| Files: | 12,203 |
| Messages: | 5,333,315 |
| Posted today: | 1 |