Re: SVD

From: sahana_paul <sahana_paul_at_nyahnyahspammersnyahnyah>
Date: Thu, 16 Feb 2006 15:52:07 +0800

Dear Madam,
                   I think now I am quite clear with the dimension, as the
confusion was due to difference between Fartran and NCL.
Still I have some cofusion. To crosscheck I calculate the A,B (time
expansion coefficients) and then calculate <A|B>, which will give diagonal
singular matrix, with singular values as the diagonal element. But
misterously I found the diagonal values are 8 time the singular values given
by the svcov_sv function. Although they should be the singular values only
and relevently ntime = 8.It is giving same error for every exapmle. I dont
find the mistake. Kindly help to find out the mistake I made. The programme
and the results are given below for your convenience.

I wait for your kind reply,

regards

sahana

begin
   ntime = 8 ;#time steps
   ncols = 4 ;#columns (stations
or grid pts) for S
   ncolz = 9 ;#columns (stations
or grid pts) for Z
   nsvd = 4 ;#svd patterns to
calculate
                                                        ;[nsvd<=min(ncols,ncolz)]

   grb_file1 = addfile("I:/ECHAM4/e4op_a2_199001.grb","r") ; read the grib
file
   grb_file2 = addfile("I:/ECHAM4/e4op_b2_199001.grb","r") ; read the grib
file

   s = grb_file1->TCC_GDS4_SFC_10(:7,:1,:1)

   z = grb_file1->2T_GDS4_SFC_10(:7,:2,:2)

   s1 = ndtooned(s)

   z1 = ndtooned(z)

   s2=onedtond(s1,(/8,4/))
   z2=onedtond(z1,(/8,9/))

   write_matrix(s2,"4F10.2",False)
   write_matrix(z2,"9F8.2",False)

   s2!0 = "time" ;name dimension for
recording

   s2!1 = "col"

   z2!0 = "time"

   z2!1 = "col"

   svLeft = new((/nsvd,ncols/),float) ; pre-allocate space
   svRight = new((/nsvd,ncolz/),float) ; reorder so time
varying fastest

   pcVar = svdcov_sv (s2(col|:,time|:),z2(col|:,time|:),nsvd,svLeft,svRight)

   print("svdcov_sv:percent variance= "+pcVar)

   x =s2(col|:,time|:)
   y = z2(col|:,time|:)

   write_matrix(svLeft,"4F10.4",False)
   write_matrix(svRight,"9F8.4",False)

    A= svLeft # x
    B= svRight # y

    write_matrix(A,"8F10.4",False)
    write_matrix(B,"8F10.4",False)

    C = (B) # transpose(A)

    write_matrix(C,"4F10.1",False)
    print("svdcov_sv:singular value= "+pcVar_at_sv)
end

The Result

Copyright (C) 1995-2004 - All Rights Reserved
 University Corporation for Atmospheric Research
 NCAR Command Language Version 4.2.0.a032
 The use of this software is governed by a License Agreement.
 See http://ngwww.ucar.edu/ncl/ for more details.
  S2:
      1.98 1.98 1.98 1.97
      0.99 0.99 0.99 0.99
      0.99 0.99 0.99 0.99
      0.99 0.99 0.87 0.91
      0.98 0.99 0.95 0.90
      0.98 0.98 0.95 0.90
      0.99 0.99 0.98 0.97
      0.99 0.99 0.99 0.97

  Z2:
  533.82 534.34 534.78 530.88 532.02 532.16 531.30 533.92 534.99
  263.86 263.93 263.97 262.14 261.84 261.42 264.82 265.93 266.01
  261.35 261.35 261.36 262.50 262.75 262.77 259.34 262.71 264.38
  258.88 258.94 258.98 259.22 259.78 261.57 250.74 255.20 258.58
  258.29 258.45 258.41 253.63 255.24 256.81 248.63 250.34 252.18
  257.84 257.86 257.89 254.99 253.26 257.68 247.89 248.08 249.33
  257.62 257.81 257.93 252.00 252.92 255.08 248.07 247.46 248.03
  256.61 256.82 256.95 251.70 252.91 254.90 250.13 250.12 252.50

(0) svdcov_sv:percent variance= 100
(1) svdcov_sv:percent variance= 7.34615e-07
(2) svdcov_sv:percent variance= 4.83665e-08
(3) svdcov_sv:percent variance= 1.8842e-11

svLeft :

   -0.5061 -0.5065 -0.4965 -0.4908
   -0.4608 -0.4918 0.3074 0.6718
    0.1490 0.1099 -0.8116 0.5541
   -0.7137 0.6997 -0.0158 0.0300

svRight :

 -0.3353 -0.3355 -0.3356 -0.3325 -0.3330 -0.3344 -0.3295 -0.3313 -0.3329
 -0.2235 -0.2136 -0.1993 -0.1653 -0.1463 -0.3949 0.5937 0.4709 0.2931
 -0.3353 -0.3506 -0.3575 0.2867 0.3087 0.2429 -0.4092 0.1644 0.4537
  0.1748 0.1987 0.1352 -0.6058 0.3996 -0.3595 -0.3869 0.1472 0.2904

A(Left Time Expansion):

   -3.9607 -1.9842 -1.9842 -1.8800 -1.9110 -1.9070 -1.9688 -1.9727
    0.0422 0.0264 0.0264 -0.0698 -0.0474 -0.0436 0.0082
0.0106
   -0.0062 0.0012 0.0012
 0.0551 -0.0138 -0.0147 -0.0054 -0.0117
   -0.0002 0.0001 0.0001 -0.0005
 0.0036 -0.0018 -0.0005 -0.0006

B(Right Time Expansion):

-1599.3895 -791.2842 -786.1675 -773.9870 -764.0369 -761.6573 -759.0292 -760.9126
    7.8449 7.5643
 3.2193 -4.1346 -5.7273 -8.0157 -6.9598 -2.4139
    1.3316 0.1815 4.5469
 4.5733 -1.8521 -2.6885 -4.9660 -2.4384
   -0.2004 0.0071 -0.4597 0.0755 1.3214 -1.7505 0.2853
0.9130

C = <A|B> = B # transpose(A) [ for NCL]

   16827.8 0.0 0.0 0.0
       0.0 1.4 0.0 0.0
       0.0 0.0 0.4 0.0
       0.0 0.0 0.0 0.0

(0) svdcov_sv:singular value= 2103.47
(1) svdcov_sv:singular value= 0.180288
(2) svdcov_sv:singular value= 0.0462604
(3) svdcov_sv:singular value= 0.000913061

----- Original Message -----
From: "Mary Haley" <haley_at_ucar.edu>
To: "sahana_paul" <sahana_paul_at_gcc.ntu.edu.tw>
Cc: "Dennis Shea" <shea_at_ucar.edu>
Sent: Thursday, February 16, 2006 12:19 PM
Subject: Re: SVD

>
> Dear Sahana,
>
> Can you elaborate on what you mean by your message below? Were
> you getting an error message from svdcov_sv, or are you saying that you
> observed a case where the wrong dimension
> size is being used?
>
> Because NCL and Fortran treat their arrays differently, from the NCL
> side the arrays are dimensioned nsvd x ncolx and nsvd x ncoly for
> svleft and svright, and ncolx x ntimes and ncoly x ntimes for athe
> input arrays.
>
> The underlying Fortran code for svdcov_sv has svleft dimensioned as
> ncolx x nsvd and svright dimensioned as ncoly x nsvd. The two input
> arrays are dimensioned ntimes x ncolx and ntimes x ncoly respectively.
>
> Perhaps this is where the confusion is coming in, with the difference
> between Fortran and NCL?
>
> --Mary
>
> On Mon, 13 Feb 2006, sahana_paul wrote:
>
>> hi,
>> While running SVD (svdcov_sv) I found the dimension of the svLeft
>> i.e. the leftsingular vector should be (nclos,nsvd); as the dimension of
>> the covariance matrix should be (nclos,ncloz), so the dimension of
>> svLeft=(ncols,nsvd)
>> sv =(nsvd,nsvd)
>> svRight=(nsvd,ncolz)
>> (ref. Bretherton et. al.
>> 1992)
>> But , I find in the pakage it is reverse dimension for svLeft. Kindly
>> solve the mistry.
>>
>> regards
>>
>> Sahana Paul.
>>
>> Postdoctoral Research Fellow
>> National Taiwan University,
>> Taiwan
>> ROC.
>
Received on Thu Feb 16 2006 - 00:52:07 MST

This archive was generated by hypermail 2.2.0 : Fri Feb 17 2006 - 15:24:09 MST