bug #44061: U and V result matrices from svds are not othogonal when rows < 9

 Submitted by: None Submitted on: Wed 21 Jan 2015 04:23:36 PM UTC Category: Octave Function Severity: 3 - Normal Priority: 5 - Normal Item Group: Incorrect Result Status: Confirmed Assigned to: None Originator Name: Svend H. SÃ¸rensen Originator Email: -email is unavailable- Open/Closed: Open Release: dev Operating System: GNU/Linux

You are not logged in

Discussion

Sun 18 Dec 2016 08:27:56 PM UTC, comment #10:

This behavior seems to be still present in Octave 4.2.0.

The script from comment #0 still gives the same results.

Hartmut <hardy>
Fri 06 Feb 2015 04:33:17 PM UTC, comment #9:

@marco: I put in another cset to make it explicit that tolower is to be used only on char inputs (http://hg.savannah.gnu.org/hgweb/octave/rev/6e2729924601).  I had done my testing on the development branch where tolower has been changed to return numeric input and only operate on char inputs.

Rik <rik5>
Thu 05 Feb 2015 09:01:09 AM UTC, comment #8:

Rik,

I think your patch should be modified, since when tmp is a number, tolower (tmp) is not defined (at least, if I try tolower(4) I see a strange symbol).

Coming back to the original bug, I have an explanation (but not yet a fix). If you call A the original matrix and try

[V,s] = eigs([sparse(4,4),A;A',sparse(4,4)],8)

you see four eigenvectors relative to 0, two of them with the last four entries 0. When you ask for 4 singular values of A, close to 4, you get the eigenvectors relative to 2 and two of the eigenvectors relative to 0, precisely the two with the last four entries 0. Try

[V,d]=eigs([sparse(4,4),A;A',sparse(4,4)],4,4)

to see it. This is precisely the command used in svds.
This is the problem, since the last entries are used to build the orthogonal matrix V. Here we get these two eigenvectors because if rows(A)<9, eigs calls eig (see call_eig in eigs.m). If eigs would call ARPACK, since it is an iterative method based on a random initial vector, you would (almost) never see such two vectors, but rather two "random" vectors spanned by the four eingenvectors relative to 0. This is a case in which a random initialization helps you. A second problem is that these two vectors, although orthogonal, are not orthognal when split into two parts. Then, you need to orthogonalize the two parts. I attached a script to demonstrate this second point.

Marco Caliari <caliari>
Thu 05 Feb 2015 02:12:50 AM UTC, comment #7:

See comment #2, I've already committed a change for the 'LA' error.  You'll need to build from a recent Mercurial tip in order to see the fix.

Rik <rik5>
Wed 04 Feb 2015 07:54:50 AM UTC, comment #6:

Andrew,

I recently upgraded to 3.8.2 and I see

V =

1
0

s =  1

which is now correct. With 3.8.1 I'm pretty sure that V was wrong, try [V,s]=eigs(rand(2),1). Anyway, the original submission and your test svds(rand(4)) fail as described.

Marco Caliari <caliari>
Wed 04 Feb 2015 07:33:43 AM UTC, comment #5:

Rik, yes Matlab does output orthonormal vectors:

U =

0.5323   -0.4654   -0.1215    0.6966
0.4654    0.5323   -0.6966   -0.1215
-0.4654   -0.5323   -0.6966   -0.1215
-0.5323    0.4654   -0.1215    0.6966

V =

-0.9978   -0.0669    0.0000    0.0000
-0.0669    0.9978    0.0000   -0.0000
-0.0000    0.0000    0.0002    1.0000
-0.0000    0.0000    1.0000   -0.0002

Marco, what is wrong with the value returned by eigs?  I get

V =

1   0

s =  1

Any (non-zero) vector is an eigenvector of the identity matrix, isn't it?  That means that the returned  V  is (the transpose of) an eigenvector corresponding to the eigenvalue returned.

However, for  n<4,  I get a failure of

[U, S, V] = svds(rand(4))
error: unrecognized value for sigma: LA
error: called from:
error:   /usr/share/octave/3.8.1/m/sparse/eigs.m at line 331, column 9
error:   /usr/share/octave/3.8.1/m/sparse/eigs.m at line 254, column 15
error:   /usr/share/octave/3.8.1/m/sparse/svds.m at line 182, column 17

in version 3.8.1

Lachlan Andrew <lachlan>
Thu 29 Jan 2015 12:53:57 PM UTC, comment #4:

1)

[U,S,V] = svds(rand(n))

gives orthogonal matrices, both Matlab and Octave (U'*U=1, not U*U').

2) But for this particular matrix, it does not work. It seems a problem of eigs when computing only some of the eigenvalues

[V,s] = eigs(eye(2),1)

V is not the matrix of eigenvectors as described in the help and as it is in Matlab. The dimensions are even wrong (and it is not just the transposed).

Marco Caliari <caliari>
Tue 27 Jan 2015 09:34:59 PM UTC, comment #3:

This may be related to interpretation, and maybe something Octave should be doing better.  If you read the documentation for svds in Octave it does not say that U or V will be orthonormal matrices such that

U*U' = 1
V*V' = 1.

The only thing it states is that

A_approx = U*S*V'

And this does in fact work with the example:

full (A)
ans =

-1  -1   0   0
-1   1   0   0
1  -1   0   0
1   1   0   0

[U, S, V] = svds (A)

U*S*V' - A
ans =

0   0   0   0
0   0   0   0
0   0   0   0
0   0   0   0

Can someone with access to Matlab run the example matrix?  According to their documentation they are supposed to output orthonormal U and V matrices.

This doesn't appear to be a problem with this particular matrix.  I tried svds on several matrices created with rand() and U and V are never orthonormal.

Rik <rik5>
Tue 27 Jan 2015 07:58:47 PM UTC, comment #2:

@Jordi: See this cset from 3 days ago on stable which fixes the uppercase argument thing (http://hg.savannah.gnu.org/hgweb/octave/rev/08dbf53f34c1).  I haven't pulled it through gui-release and on to dev yet.  I figured there would be something else that would flush that cset down.

Rik <rik5>
Mon 26 Jan 2015 04:43:42 PM UTC, comment #1:

Something is definitely fishy here... if nothing else, I see a failure with svds (A), because eigs () does not accept uppercase arguments.

Wed 21 Jan 2015 04:23:36 PM UTC, original submission:

octave:2> A = spconvert(two_vars_neg_in)
A =

Compressed Column Sparse (rows = 4, cols = 4, nnz = 8 [50%])

(1, 1) -> -1
(2, 1) -> -1
(3, 1) ->  1
(4, 1) ->  1
(1, 2) -> -1
(2, 2) ->  1
(3, 2) -> -1
(4, 2) ->  1

octave:3>  [U,S,V]= svds(A,rows(A),rows(A))
U =

-0.50000   0.50000   1.00000   0.00000
0.50000   0.50000   0.00000   1.00000
-0.50000  -0.50000   0.00000   1.00000
0.50000  -0.50000   1.00000   0.00000

S =

Diagonal Matrix

2   0   0   0
0   2   0   0
0   0   0   0
0   0   0   0

V =

-0.00000  -1.00000   0.00000   0.00000
1.00000   0.00000   0.00000   0.00000
0.00000   0.00000   0.00000   0.00000
0.00000   0.00000   0.00000   0.00000

octave:4> #[U,S,V] = svd (full(A) )
octave:4> ss = diag(S)
ss =

2
2
0
0

octave:5> nnull=sum(ss(:)==0)
nnull =  2
octave:6> V'*V  ## this should result en the identity matrix but does not !!
ans =

1.00000   0.00000   0.00000   0.00000
0.00000   1.00000   0.00000   0.00000
0.00000   0.00000   0.00000   0.00000
0.00000   0.00000   0.00000   0.00000

octave:7> U'*U  ## here again we expect to see the identity matrix
ans =

1.00000   0.00000   0.00000   0.00000
0.00000   1.00000   0.00000   0.00000
0.00000   0.00000   2.00000   0.00000
0.00000   0.00000   0.00000   2.00000

Anonymous

Attached Files

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Comment:

Attached Files
file #32999:  script.m added by caliari (370B - text/x-objcsrc)
file #32889:  two_vars_neg_in.sparse added by None (121B - application/octet-stream)
file #32890:  svd-2-var.oct added by None (178B - application/octet-stream)

Dependencies

Depends on the following items: None found

Items that depend on this one: None found

Carbon-Copy List
• -email is unavailable- added by hardy (Posted a comment)
• -email is unavailable- added by lachlan (Posted a comment)
• -email is unavailable- added by caliari (Posted a comment)
• -email is unavailable- added by rik5 (Updated the item)
• -email is unavailable- added by jordigh (Posted a comment)
• -email is unavailable- added by None (Submitted the item)
•

Do you think this task is very important?
This task has 0 encouragements so far.

Only project members can vote.