Opened 7 years ago

Closed 7 years ago

#13140 closed defect (fixed)

OS X Lion doctest failures for double dense QR decomposition

Reported by: jdemeyer Owned by: jason, was
Priority: major Milestone: sage-5.2
Component: linear algebra Keywords:
Cc: rbeezer Merged in: sage-5.2.beta1
Authors: John Palmieri, Rob Beezer Reviewers: Rob Beezer, John Palmieri
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: #10795 Stopgaps:

Description (last modified by jhpalmieri)

The following doctest failure on OS X Lion is introduced by #10795:

sage -t  "devel/sage/sage/matrix/matrix_double_dense.pyx"  
**********************************************************************
File "/Users/palmieri/Desktop/Sage_stuff/sage_builds/clean/sage-5.1.beta5/devel/sage/sage/matrix/matrix_double_dense.pyx", line 2477:
    sage: Q.round(6).zero_at(10^-6)
Expected:
    [-0.458831  0.126051 -0.381212 -0.394574  -0.68744]
    [-0.458831  -0.47269  0.051983  0.717294 -0.220963]
    [ 0.229416 -0.661766 -0.661923 -0.180872  0.196411]
    [ 0.688247 -0.189076  0.204468   0.09663 -0.662889]
    [-0.229416 -0.535715  0.609939 -0.536422  0.024551]
Got:
    [-0.458831 -0.126051  0.381212 -0.394574  -0.68744]
    [-0.458831   0.47269 -0.051983  0.717294 -0.220963]
    [ 0.229416  0.661766  0.661923 -0.180872  0.196411]
    [ 0.688247  0.189076 -0.204468   0.09663 -0.662889]
    [-0.229416  0.535715 -0.609939 -0.536422  0.024551]
**********************************************************************
File "/Users/palmieri/Desktop/Sage_stuff/sage_builds/clean/sage-5.1.beta5/devel/sage/sage/matrix/matrix_double_dense.pyx", line 2483:
    sage: R.round(6).zero_at(10^-6)
Expected:
    [ 4.358899 -0.458831 13.076697  6.194225  2.982405]
    [      0.0 -1.670172 -0.598741   1.29202 -6.207997]
    [      0.0       0.0 -5.444402 -5.468661  0.682716]
    [      0.0       0.0       0.0  1.027626   -3.6193]
    [      0.0       0.0       0.0       0.0  0.024551]
Got:
    [ 4.358899 -0.458831 13.076697  6.194225  2.982405]
    [      0.0  1.670172  0.598741  -1.29202  6.207997]
    [      0.0       0.0  5.444402  5.468661 -0.682716]
    [      0.0       0.0       0.0  1.027626   -3.6193]
    [      0.0       0.0       0.0       0.0  0.024551]
**********************************************************************
File "/Users/palmieri/Desktop/Sage_stuff/sage_builds/clean/sage-5.1.beta5/devel/sage/sage/matrix/matrix_double_dense.pyx", line 2511:
    sage: Q.round(6).zero_at(10^-6)
Expected:
    [             -0.730297  0.207057 + 0.538347*I -0.246305 + 0.076446*I   0.238162 - 0.10366*I]
    [              0.091287 -0.207057 - 0.377878*I -0.378656 + 0.195222*I  0.701244 - 0.364371*I]
    [  0.63901 + 0.091287*I  0.170822 + 0.667758*I  0.034115 - 0.040902*I  0.314017 - 0.082519*I]
    [ 0.182574 + 0.091287*I  -0.036235 + 0.07247*I -0.863228 - 0.063228*I -0.449969 - 0.011612*I]
Got:
    [             -0.730297  0.207057 + 0.538347*I  0.246305 - 0.076446*I   0.238162 - 0.10366*I]
    [              0.091287 -0.207057 - 0.377878*I  0.378656 - 0.195222*I  0.701244 - 0.364371*I]
    [  0.63901 + 0.091287*I  0.170822 + 0.667758*I -0.034115 + 0.040902*I  0.314017 - 0.082519*I]
    [ 0.182574 + 0.091287*I  -0.036235 + 0.07247*I  0.863228 + 0.063228*I -0.449969 - 0.011612*I]
**********************************************************************
File "/Users/palmieri/Desktop/Sage_stuff/sage_builds/clean/sage-5.1.beta5/devel/sage/sage/matrix/matrix_double_dense.pyx", line 2516:
    sage: R.round(6).zero_at(10^-6)
Expected:
    [             10.954451            -1.917029*I   5.385938 - 2.19089*I  -0.273861 - 2.19089*I]
    [                   0.0               4.829596 -0.869638 - 5.864879*I  0.993872 - 0.305409*I]
    [                   0.0                    0.0             -12.001608  0.270953 - 0.442063*I]
    [                   0.0                    0.0                    0.0               1.942964]
Got:
    [             10.954451            -1.917029*I   5.385938 - 2.19089*I  -0.273861 - 2.19089*I]
    [                   0.0               4.829596 -0.869638 - 5.864879*I  0.993872 - 0.305409*I]
    [                   0.0                    0.0              12.001608 -0.270953 + 0.442063*I]
    [                   0.0                    0.0                    0.0               1.942964]
**********************************************************************
File "/Users/palmieri/Desktop/Sage_stuff/sage_builds/clean/sage-5.1.beta5/devel/sage/sage/matrix/matrix_double_dense.pyx", line 2548:
    sage: R.round(6).zero_at(10^-6)
Expected:
    [-5.567764   2.69408  -2.69408]
    [      0.0 -3.569585  3.569585]
    [      0.0       0.0       0.0]
    [      0.0       0.0       0.0]
Got:
    [ 5.567764  -2.69408   2.69408]
    [      0.0  3.569585 -3.569585]
    [      0.0       0.0       0.0]
    [      0.0       0.0       0.0]
**********************************************************************

Apply trac_13140_double_dense.patch and trac_13140-rescale-inplace.patch.

Attachments (2)

trac_13140_double_dense.patch (7.5 KB) - added by jhpalmieri 7 years ago.
trac_13140-rescale-inplace.patch (673 bytes) - added by rbeezer 7 years ago.

Download all attachments as: .zip

Change History (11)

comment:1 Changed 7 years ago by rbeezer

  • Cc rbeezer added

Appears to be eigenvectors that differ by a sign. Normalizing to a leading positive entry will make them unique and predictable.

comment:2 Changed 7 years ago by jdemeyer

  • Milestone changed from sage-5.1 to sage-5.2
  • Priority changed from blocker to major

comment:3 follow-up: Changed 7 years ago by jhpalmieri

  • Authors set to John Palmieri
  • Status changed from new to needs_review

Here's a possible patch. Rob, what do you think?

comment:4 in reply to: ↑ 3 ; follow-up: Changed 7 years ago by rbeezer

Replying to jhpalmieri:

Here's a possible patch. Rob, what do you think?

Dear John,

Thanks for working on this. This should be my mess to cleanup, but I have not had a chance to get to it.

Some comments. I get 7 doctest failures. This is on 5.1.beta3, but the offending #10795 applied cleanly, along with your #13140.

It looks to me like the following stanza in _normalize_columns() finds the first negative entry of the column, and adjusts that. Which would not be the leading (non-zero) entry.

for j from 0 <= j < M.ncols():
    for i from 0 <= i < M.column(j).degree():
        if M.column(j)[i] < 0:
            M = M.with_rescaled_col(j, -1)
            break

If I replace this with

for j from 0 <= j < M.ncols():
    for i from 0 <= i < M.column(j).degree():
        if M.column(j)[i] != 0:
            if M.column(j)[i] < 0:
                M = M.with_rescaled_col(j, -1)
            break

then the number of failures goes down to 4. Notice this does not carefully consider the possibility of -0.0 before the desired leading entry (can we just test on that literal?).

Line 2517 (new line number according to patch) has a .normalize_columns() without an underscore.

With my change above I think any remaining problems (barring one) are just +/- 0 mixups. It might make sense to normalize first, then round/clip to six places.

Most critically - Q*R should be the original matrix. So if we adjust a column of Q, then we need to adjust a row of R before testing Q*R == A. We could test this property before normalizing. Then a companion _normalize_rows() might be indicated before displaying Q and R. But still, the rows of R to normalize may not be the same as the columns of Q to normalize.

I'll do my best to stick with this one, but the next three days are full of family obligations (taking care of the parental generation!).

Thanks again for working on this. It'd be good to have a pair of these normalization routines available.

Rob

Changed 7 years ago by jhpalmieri

comment:5 in reply to: ↑ 4 ; follow-up: Changed 7 years ago by jhpalmieri

Hi Rob,

Here's another attempt.

Replying to rbeezer:

Some comments. I get 7 doctest failures. This is on 5.1.beta3, but the offending #10795 applied cleanly, along with your #13140.

When applied to 5.1.beta6, the new patch passes tests for me on OS X Lion, sage.math, and hawk (OpenSolaris).

It looks to me like the following stanza in _normalize_columns() finds the first negative entry of the column, and adjusts that. Which would not be the leading (non-zero) entry.

for j from 0 <= j < M.ncols():
    for i from 0 <= i < M.column(j).degree():
        if M.column(j)[i] < 0:
            M = M.with_rescaled_col(j, -1)
            break

You're right, thanks, I've fixed it.

then the number of failures goes down to 4. Notice this does not carefully consider the possibility of -0.0 before the desired leading entry (can we just test on that literal?).

Within Sage, RDF(-0.0) prints as "-0.0", so I think it's the right thing to look at. RDF(-0.0) == 0 returns True while RDF(-0.0) < 0 returns False, which suggests that we're okay, but I'm not absolutely sure.

Line 2517 (new line number according to patch) has a .normalize_columns() without an underscore.

Fixed.

With my change above I think any remaining problems (barring one) are just +/- 0 mixups. It might make sense to normalize first, then round/clip to six places.

Right, good idea.

Most critically - Q*R should be the original matrix. So if we adjust a column of Q, then we need to adjust a row of R before testing Q*R == A. We could test this property before normalizing. Then a companion _normalize_rows() might be indicated before displaying Q and R. But still, the rows of R to normalize may not be the same as the columns of Q to normalize.

In the new patch, when printing Q or R for doctests, we normalize them. When multiplying them, we don't normalize them. Then we don't have to worry about compatibility of normalizations. I'm also now normalizing the rows of R and the columns of Q, since those seem like the most natural choices.

comment:6 in reply to: ↑ 5 Changed 7 years ago by rbeezer

Replying to jhpalmieri:

Here's another attempt.

Sounds good and looks good. It'll be Tuesday probably before I can give it a proper review, but I'll get to it as quickly as I can then.

Rob

Changed 7 years ago by rbeezer

comment:7 Changed 7 years ago by rbeezer

  • Dependencies set to #10795
  • Reviewers set to Rob Beezer

Everything checks out fine on version 5.1.rc1 and all looks good. Thanks for putting this together, maybe QR decomposition will finally be in good order.

One final idea. I made a slightly different version of this routine ("inplace"), using the matrix routine that scales the column without making a copy of the matrix. Runs about 10% faster on a moderately large matrix.

sage: A = random_matrix(RDF, 1000)
sage: timeit("A._normalize_columns()")
5 loops, best of 3: 5.38 s per loop
sage: timeit("A._normalize_columns_inplace()")
5 loops, best of 3: 4.84 s per loop

Patch implements the suggested change, and the relevant tests pass when applied. If you like it, then go ahead and add it to the list of patches and flip this all to positive review. Thanks again.

comment:8 Changed 7 years ago by jhpalmieri

  • Authors changed from John Palmieri to John Palmieri, Rob Beezer
  • Description modified (diff)
  • Reviewers changed from Rob Beezer to Rob Beezer, John Palmieri
  • Status changed from needs_review to positive_review

Looks good to me.

comment:9 Changed 7 years ago by jdemeyer

  • Merged in set to sage-5.2.beta1
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.