Re: Behaviour of svd_lapack

Lunde, Bruce N CIV NAVOCEANO, NP1 wrote:

> I have some questions about the behaviour
> or implementation of the NCL routine
>
> svd_lapack
>

You say you have questions but do not ask any here but rather make a
series of observations. I can only say whether your statements are true
or not.

> I have been testing it against the SVD
> results of MATLAB and Octave (a free MATLAB
> clone). The results of MATLAB and Octave
> agree (not surprising), but NCL is not
> "quite" the same.
>
> Bottom line is that svd_lapack
> (0) needs the inputs (U,V) to be float or double,
>

True.

> (1) seems to not like nrow.GT.ncol ,
>

True, and I do not see any reason for this. The underlying LAPACK
routines can calculate singular values and singular vectors (left and
right) for any M,N > 0.

> (2) needs the input matrix transposed to get
> "close to" the same results that MATLAB does,
>

This is purely a notational issue. Ordinarily we talk about MxN matrices
which have M rows and N columns. These are referred to in NCL as NxM
arrays. Yeah, I agree, it's quirky if you come from a Fortran / MATLAB
background.

> (3) does not return the full U matrix,
>

Only the first 2 columns of U have meaning, even in the LAPACK and
MATLAB versions of this problem. These are the left singular vectors. So
in this sense, NCL *does* return the full U matrix.

When calling the underlying LAPACK routines, there are two
single-character arguments called JOBU and JOBV. The answer you got with
MATLAB sets these two arguments to "A" while NCL sets them to "S" which
means that only the first min(M,N) columns of U are returned. (NCL does
not implement the "A" option through the user interface.)

> (4) returns a transposed U matrix but correct V
> matrix when compared to MATLAB (after
> tranposing the input matrix),

Again, notation only.

> (5) but the numbers it does return DO MATCH those of
> MATLAB's svd.
>

Well, they'd better match or something is seriously wrong.

> a = (/ (/1,2/), (/3,4/), (/5,6/), (/7,8/) /)*1.0d0
>
> ; print(dimsizes(a)) gives 4 2 ... matching MATLAB
>

But you need 2 4 to match MATLAB.

My version of the NCL program (which, note, does get the right answer
when you interpret rows and columns correctly):

  a = new((/2,4/),"double")
  a(0,0) = 1.0d0
  a(1,0) = 2.0d0
  a(0,1) = 3.0d0
  a(1,1) = 4.0d0
  a(0,2) = 5.0d0
  a(1,2) = 6.0d0
  a(0,3) = 7.0d0
  a(1,3) = 8.0d0

  print(dimsizes(a))

  dima = dimsizes(a)
  nca = dima(0)
  nra = dima(1)

  nrcu = min( (/nra,nca/) ) ; 4, option S

  u = new ( (/nrcu,nra /) , typeof(a) ) ; pre-allocate space
  v = new ( (/nrcu,nrcu/) , typeof(a) ) ; for returned arrays

  s = svd_lapack (a, "S" , "S", 0, u, v)

More questions, comments or observations are welcome,

Andy
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Jul 02 2009 - 14:10:19 MDT