[6164d0]: mp_mysvd.m Maximize Restore History

Download this file

mp_mysvd.m    297 lines (293 with data), 7.4 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
function [U,S,V]=mp_mysvd(A,dummy)
%function [U,S,V]=mysvd(A,dummy)
%svd: singular value decomposition
%Outputs should verify that U*S*V'==A
%The case with two arguments is considered as an "economy size" request
% Based upon FORTRAN 77 routine at NUMERICAL RECIPES, monkey-modified to run in Matlab
% also taking into consideration the C implementation (due to the GOTO hell!)
% Given a matrix A(1:m,1:n),with physical dimensions mp by np this routine computes its
% singular value decomposition,A = U * W * V' .The matrix U replaces A on output.The
% diagonal matrix of singular values W is output as a vector w(1:n). The matrix V (not the
% transpose V') is output as v(1:n,1:n).
% NOTICE that direct comparison with the matlab svd is possible only for the singular values,
% but for the left and right matrix you should be prepared to consider some possible
% sign differences. The verification should be instead that norm(U*S*V'-A) is small, and
% that norm(S-Smatlab) is also small
%
%source in FORTRAN 77 at http://cc.oulu.fi/~tf/tiedostot/pub/nrf/svdcmp.f
%source in C at http://www.cs.utah.edu/~tch/classes/cs4550/code/Numerical-Recipes/sources-used/svdcmp.c
[m,n]=size(A);
a_zero=A(1)*0;%a zero of the same data type as input matrix A
a_one=1+a_zero;%a one of the same data type as input matrix A
g=a_zero;
scale=a_zero;
anorm=a_zero;
for i=1:n %do 25 i=1,n
L=i+1;
rv1(i)=scale*g;
g=a_zero;
s=a_zero;
scale=a_zero;
if(i<=m)
for k=i:m %do 11 k=i,m
scale=scale+abs(A(k,i));
end %enddo 11
if (scale~=0)
for k=i:m%do 12 k=i,m
A(k,i)=A(k,i)/scale;
s=s+A(k,i)*A(k,i);
end %enddo 12
f=A(i,i);
g=-signtransf(sqrt(s),f);
h=f*g-s;
A(i,i)=f-g;
if i~=n
for j=L:n %do 15 j=l,n
s=a_zero;
for k=i:m %do 13 k=i,m
s=s+A(k,i)*A(k,j);
end %enddo 13
f=s/h;
for k=i:m %do 14 k=i,m
A(k,j)=A(k,j)+f*A(k,i);
end %enddo 14
end %enddo 15
end
for k=i:m %do 16 k=i,m
A(k,i)=scale*A(k,i);
end %enddo 16
end
end
w(i)=scale *g;
g=a_zero;
s=a_zero;
scale=a_zero;
if((i<=m) & (i~=n))
for k=L:n %do 17 k=l,n
scale=scale+abs(A(i,k));
end %enddo 17
if(scale~=a_zero)
for k=L:n %do 18 k=l,n
A(i,k)=A(i,k)/scale;
s=s+A(i,k)*A(i,k);
end %enddo 18
f=A(i,L);
g=-signtransf(sqrt(s),f);
h=f*g-s;
A(i,L)=f-g;
for k=L:n %do 19 k=l,n
rv1(k)=A(i,k)/h;
end %enddo 19
for j=L:m %do 23 j=l,m
s=a_zero;
for k=L:n %do 21 k=l,n
s=s+A(j,k)*A(i,k);
end %enddo 21
for k=L:n %do 22 k=l,n
A(j,k)=A(j,k)+s*rv1(k);
end %enddo 22
end %enddo 23
for k=L:n %do 24 k=l,n
A(i,k)=scale*A(i,k);
end %enddo 24
end
end
anorm=max(anorm,(abs(w(i))+abs(rv1(i))));
end %enddo 25
v=repmat(a_zero,n,n);%Preallocate space
for i=n:-1:1 %do 32 i=n,1,-1 Accumulation of right-hand transformations.
if(i<n)
if(g~=a_zero)
for j=L:n %do 26 j=l,n Double division to avoid possible underflow.
v(j,i)=(A(i,j)/A(i,L))/g;
end %enddo 26
for j=L:n %do 29 j=l,n
s=a_zero;
for k=L:n %do 27 k=l,n
s=s+A(i,k)*v(k,j);
end %enddo 27
for k=L:n %do 28 k=l,n
v(k,j)=v(k,j)+s*v(k,i);
end %enddo 28
end %enddo 29
end
for j=L:n %do 31 j=l,n
v(i,j)=a_zero;
v(j,i)=a_zero;
end %enddo 31
end
v(i,i)=1+a_zero;
g=rv1(i);
L=i;
end %enddo 32
for i=min(m,n):-1:1 %do 39 i=min(m,n),1,-1 Accumulation of left-hand transformations.
L=i+1;
g=w(i);
for j=L:n %do 33 j=l,n
A(i,j)=a_zero;
end %enddo 33
if(g~=a_zero)
g=a_one/g;
for j=L:n %do 36 j=l,n
s=a_zero;
for k=L:m %do 34 k=l,m
s=s+A(k,i)*A(k,j);
end %enddo 34
f=(s/A(i,i))*g;
for k=i:m %do 35 k=i,m
A(k,j)=A(k,j)+f*A(k,i);
end %enddo 35
end %enddo 36
for j=i:m %do 37 j=i,m
A(j,i)=A(j,i)*g;
end %enddo 37
else
for j=i:m %do 38 j= i,m
A(j,i)=a_zero;
end %enddo 38
end
A(i,i)=A(i,i)+a_one;
end %enddo 39
for k=n:-1:1 %do 49 k=n,1,-1
%Diagonalization of the bidiagonal form:Loop over
%singular values,and over allowed iterations.
for its=1:30 %do 48 its=1,30
flag=1;%default action
for L=k:-1:1 %do 41 l=k,1,-1 Test for splitting.
nm=L-1; %Note that rv1(1) is always zero.
if((abs(rv1(L))+anorm)==anorm)
flag=0; %goto 2
break
elseif((abs(w(nm))+anorm)==anorm)
break
end
end %enddo 41
if flag==1
%1 label 1
c=a_zero; % Cancellation of rv1(l),if l > 1
s=a_one;
for i=L:k %do 43 i=l,k
f=s*rv1(i);
rv1(i)=c*rv1(i);
if((abs(f)+anorm)==anorm)
break
end
g=w(i);
h=pythag(f,g);
w(i)=h;
h=a_one/h;
c= (g*h);
s=-(f*h);
for j=1:m %do 42 j=1,m
y=A(j,nm);
z=A(j,i);
A(j,nm)=(y*c)+(z*s);
A(j,i)=-(y*s)+(z*c);
end %enddo 42
end %enddo 43
end %del goto 1
%2
z=w(k);
if(L==k)%then Convergence.
if(z<=a_zero)% Singular value is made nonnegative.
w(k)=-z;
for j=1:n %do 44 j=1,n
v(j,k)=-v(j,k);
end %enddo 44
end
break
end
if(its==30)
error('no convergence in 30 iterations of svdcmp')
end
x=w(L);% Shift from bottom 2-by-2 minor.
nm=k-1;
y=w(nm);
g=rv1(nm);
h=rv1(k);
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g=pythag(f,a_one);
f=((x-z)*(x+z)+h*((y/(f+signtransf(g,f)))-h))/x;
c=a_one; %Next QR transformation:
s=a_one;
for j=L:nm %do 47 j=l,nm
i=j+1;
g=rv1(i);
y=w(i);
h=s*g;
g=c*g;
z=pythag(f,h);
rv1(j)=z;
c=f/z;
s=h/z;
f= (x*c)+(g*s);
g=-(x*s)+(g*c);
h=y*s;
y=y*c;
for jj=1:n %do 45 jj=1,n
x=v(jj,j);
z=v(jj,i);
v(jj,j)= (x*c)+(z*s);
v(jj,i)=-(x*s)+(z*c);
end %enddo 45
z=pythag(f,h);
w(j)=z;% Rotation can be arbitrary if z =0
if(z~=a_zero)
z=a_one/z;
c=f*z;
s=h*z;
end
f= (c*g)+(s*y);
x=-(s*g)+(c*y);
for jj=1:m %do 46 jj=1,m
y=A(jj,j);
z=A(jj,i);
A(jj,j)= (y*c)+(z*s);
A(jj,i)=-(y*s)+(z*c);
end %enddo 46
end %enddo 47
rv1(L)=a_zero;
rv1(k)=f;
w(k)=x;
end %enddo 48
%3 continue
end %enddo 49
[dummy,I]=sort(-w(1:n));%Sort in decreasing order
S=diag(w(I));
if nargin<2
%fill in with zeros in order to satisfy dimension requirements
if n<size(A,1),S(size(A,1),n)=a_zero;end
if n<size(A,2),S(n,size(A,2))=a_zero;end
%Create a square U matrix
dummy=A(:,I);
[U,dummy]=qr(dummy);
%correct some signs
q=find(diag(dummy)<0);U(:,q)=-U(:,q);
else
%go on with the economy size version
U=A(:,I);
end
V=v(1:n,I);
if nargout<=1,U=w(I);U=U(:);end
return
function py=pythag(a,b)
%REAL a,b,pythag
%Computes (a^2 + b^2 )^1/2 without destructive under flow or over flow.
absa=abs(a);
absb=abs(b);
if(absa>absb)
py=absa*sqrt(1.+(absb/absa)^2);
else
if(absb==0)
py=absb;
else
py=absb*sqrt(1.+(absa/absb)^2);
end
end
function c=signtransf(A1,A2)
%C = SIGNTRANSF(A1,A2) performs sign transfer: if A2 is negative the result is -abs(A1),
%if A2 is zero or positive the result is abs(A1).
%It is the fortran equivalent of SIGN(A1,A2)
% (see http://www.math.hawaii.edu/lab/197/fortran/fort4.htm#sign)
c=abs(A1);
q=find(A2<0);c(q)=-abs(A1(q));