Hi All,
Can anyone tell me why the singular value decomposition (svd_lapack) below
yields different eigenvalues than the eigenanalysis (eofunc, eofunc_ts),
while both output the exact same EOFs and PCs? Can't see what I'm doing
wrong. Code and output are below.
Thanks!
Kerrie
;######################################################################################
;### data array A=(/ntime,nspace/) ###
A=(/ (/-2.,1./),(/7.,-8./),(/-3.,0./),(/5.,-9./) /)
A!0="time"
A!1="space"
ntime=4 ;### number of rows, time dimension 0 ###
nspace=2 ;### number of columns, space dimension 1 ###
Aprime = dim_rmvmean_n(A,0)
;### singular value decomposition of Aprime ###
;### Aprime=(u)(s)(vtranspose)-----> (4x2)=(4x2)(2x2)(2x2) ###
u = new ( (/ntime,nspace/) , float ) ;### columns
of u are PC eigenvectors ###
vtranspose = new ( (/nspace,nspace/) , float ) ;###
rows of vtranspose are EOF eigenvectors ###
s = svd_lapack (Aprime, "S" , "S", 0, vtranspose, u) ;###
s^2=eigenvalues
;### s is returned with only the non-zero (diagonal) values
print("for SVD")
print("U:")
write_matrix (u, "2f10.4", False)
print("V:")
write_matrix (vtranspose, "2f10.4", False)
print("lambda:")
print(s^2)
At=transpose(A) ;### make time the rightmost dimension ###
eof=eofunc(At,2,False) ;### rows of eof are EOFs ###
eoftranspose=transpose(eof)
pc=eofunc_ts_Wrap(At,eof,False) ;### rows of pc are PCs ###
pctranspose=transpose(pc) ;### for similarity to SVD matrices,
now columns of pc are PCs ###
pc_colsum=dim_sum_n(pctranspose^2,0) ;### square each component
and sum the columns ###
pctranspose(:,0)=pctranspose(:,0)/sqrt(pc_colsum(0)) ;###
normalize PC1
pctranspose(:,1)=pctranspose(:,1)/sqrt(pc_colsum(1)) ;###
normalize PC2
print("for eofunc")
print("U:")
write_matrix (pctranspose, "2f10.4", False)
print("V:")
write_matrix (eof, "2f10.4", False)
print("lambda:")
print(eof@eval)
;######################################################################################
OUPUT IS:
(0) for SVD
(0) U:
0.5010 0.4050
-0.5261 0.5745
0.4983 -0.3748
-0.4732 -0.6046
(0) V:
-0.6898 0.7240
0.7240 0.6898
(0) lambda:
(0) 153.4626
(1) 3.287447
(0)
(0)
(0) for eofunc
(0) U:
0.5010 0.4050
-0.5261 0.5745
0.4983 -0.3748
-0.4732 -0.6046
(0) V:
-0.6898 0.7240
0.7240 0.6898
(0) lambda:
(0) 51.15419
(1) 1.095816
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Nov 29 14:14:54 2011
This archive was generated by hypermail 2.1.8 : Wed Nov 30 2011 - 19:52:47 MST