Opened 10 years ago

Closed 10 years ago

#10802 closed enhancement (fixed)

Singular values of matrices over CDF

Reported by: rbeezer Owned by: jason, was
Priority: minor Milestone: sage-4.7.2
Component: linear algebra Keywords: sd32
Cc: Merged in: sage-4.7.2.alpha4
Authors: Rob Beezer Reviewers: Martin Raum, Jeroen Demeyer
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: #10837 Stopgaps:

Status badges

Description (last modified by rbeezer)

This method serves two purposes:

(a) a convenience for getting singular values and clipping small values to zero.

(b) a starting point for reasonable quality control on numerical approximations from CDF matrices. This will make a rank computation easy (though not fool-proof) and this can be used in other situations as a consistency check. See #7392.


Depends:

  1. #10837

Apply:

  1. trac_10802-singular-values-CDF-v3.patch
  2. trac-10802-singular-values-CDF-review-v2.patch
  3. trac_10802-doctest-noise.patch

Attachments (6)

trac_10802-singular-values-CDF.patch (8.1 KB) - added by rbeezer 10 years ago.
trac_10802-singular-values-CDF-v2.patch (10.2 KB) - added by rbeezer 10 years ago.
trac-10802-singular-values-CDF-review.patch (902 bytes) - added by mraum 10 years ago.
trac-10802-singular-values-CDF-review-v2.patch (1.1 KB) - added by rbeezer 10 years ago.
trac_10802-singular-values-CDF-v3.patch (10.2 KB) - added by rbeezer 10 years ago.
Rebase
trac_10802-doctest-noise.patch (1.3 KB) - added by rbeezer 10 years ago.
Numerical noise in doctests

Download all attachments as: .zip

Change History (30)

Changed 10 years ago by rbeezer

comment:1 Changed 10 years ago by rbeezer

  • Authors set to Rob Beezer
  • Status changed from new to needs_review

comment:2 Changed 10 years ago by mraum

  • Description modified (diff)
  • Reviewers set to Martin Raum
  • Status changed from needs_review to positive_review

This looks very good to me; No complains at all. All test pass.

comment:3 Changed 10 years ago by jdemeyer

  • Status changed from positive_review to needs_work

This patch conflicts with #10837.

comment:4 follow-up: Changed 10 years ago by mraum

Robert, I think it is a good idea to rebase this patch to #10837. It's just very many changes at the same time, so this had to happen.

But also mind that the patch bot complains. None of the issues showed up on my computer, but it seems like the list tests are obviously wrong. Also the numerical noise needs to be addressed. You might want to compare against a bound 10-15 or so.

comment:5 in reply to: ↑ 4 Changed 10 years ago by jdemeyer

Replying to mraum:

But also mind that the patch bot complains. None of the issues showed up on my computer, but it seems like the list tests are obviously wrong. Also the numerical noise needs to be addressed. You might want to compare against a bound 10-15 or so.

Indeed, this also need to be fixed.

comment:6 Changed 10 years ago by rbeezer

Thanks, Martin and Jeroen, for the help and suggestions. I'm building a 4.7.alpha3 right now (and working on another project) and will get to this (and other conflicts) in a few hours from now.

Rob

Changed 10 years ago by rbeezer

comment:7 Changed 10 years ago by rbeezer

  • Description modified (diff)
  • Status changed from needs_work to needs_review

v2 patch rebases to depend on #10837 on 4.7.alpha3 and reworks the doctests.

One new doctest at the beginning that is a full-rank matrix (no zero singular values!) and should be stable numerically. includes a confidence building test on correctness.

All the other matrices stay the same, including the difficult Hilbert matrix. Using rounding liberally, plus ellipsis (...) and inequalities in a few places to get everything to behave, without obscuring the real nature of the method.

No changes to the code. New doctests pass on my Intel i7-2600 machine and on sage.math.

The changes to the doctests should be reviewed.

comment:8 Changed 10 years ago by chapoton

  • Dependencies set to #10837

comment:9 Changed 10 years ago by mraum

  • Description modified (diff)

I propose only one tiny change to the doctest involving the condition: The exact value is

sage: c2 = sqrt(sum(map(operator.mul, B.inverse().list(), B.inverse().list())))
sage: c1 = sqrt(sum(map(operator.mul, B.list(), B.list())))
sage: c1/c2
1/338581401811897828348589105304331884576*sqrt(805082999949990227/30)*sqrt(158090646384497588170245117213)
sage: N(c1*c2)
1.75177509695064e16

so indeed only the first two digits of what you computed are right. I changed the doctest accordingly.

If also modified one if clause; it now checks eps is None.

If you approve the changes, we can give this a "virtual positive review", depending on the approval of #10837.

comment:10 Changed 10 years ago by rbeezer

Thanks for the improved accuracy on that test. Looks good.

Yes, let's call this "virtual positive review" and I will wait on #10837 before actually flipping the ticket.

Thanks for your careful work on these. The improvements are welcome.

Rob

comment:11 Changed 10 years ago by rbeezer

  • Description modified (diff)
  • Status changed from needs_review to positive_review

Review patch looks good, and passes long tests, and #10837 is ready-to-go. So "positive review" here.

I made a v2 version of the review patch with a commit string, and Martin as the user.

comment:12 Changed 10 years ago by jdemeyer

  • Milestone changed from sage-4.7.1 to sage-4.7.2

comment:13 Changed 10 years ago by jdemeyer

  • Merged in set to sage-4.7.2.alpha2
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:14 Changed 10 years ago by jdemeyer

  • Resolution fixed deleted
  • Status changed from closed to new

Sorry, too quick to close this ticket. One of the patches applies only with fuzz 2 (probably because of #10839), so this should be rebased to sage-4.7.2.alpha1 + #10837.

comment:15 Changed 10 years ago by jdemeyer

  • Status changed from new to needs_review

comment:16 Changed 10 years ago by jdemeyer

  • Status changed from needs_review to needs_work

comment:17 Changed 10 years ago by jdemeyer

  • Merged in sage-4.7.2.alpha2 deleted

Changed 10 years ago by rbeezer

Rebase

Changed 10 years ago by rbeezer

Numerical noise in doctests

comment:18 follow-up: Changed 10 years ago by rbeezer

  • Description modified (diff)
  • Status changed from needs_work to needs_review

4.7.2.alpha1 + #10837

"v3" is a very trivial rebase for fuzzieness.

Two doctests now fail (did ATLAS change?). "doctest-noise" patch brutally makes these failures go away, but will necessitate a review.

comment:19 in reply to: ↑ 18 ; follow-up: Changed 10 years ago by jdemeyer

Replying to rbeezer:

Two doctests now fail (did ATLAS change?). "doctest-noise" patch brutally makes these failures go away, but will necessitate a review.

On which system did you see this? On Linux x86_64 I did not get such a failure. It is true that ATLAS did change in sage-4.7.2.alpha0.

comment:20 in reply to: ↑ 19 Changed 10 years ago by rbeezer

Replying to jdemeyer:

On which system did you see this? On Linux x86_64 I did not get such a failure. It is true that ATLAS did change in sage-4.7.2.alpha0.

Jeroen,

Here are the failures (these are the only ones on sage/matrix/matrix_double_dense.pyx).

rob@pearl:/sage/sage-4.7.2.alpha1/devel/sage$ ../../sage -t sage/matrix/matrix_double_dense.pyx
sage -t  "devel/sage-main/sage/matrix/matrix_double_dense.pyx"
**********************************************************************
File "/sage/sage-4.7.2.alpha1/devel/sage-main/sage/matrix/matrix_double_dense.pyx", line 968:
    sage: sv[2:3]
Expected:
    [2.92724029018e-16]
Got:
    [2.01161346159e-16]
**********************************************************************
File "/sage/sage-4.7.2.alpha1/devel/sage-main/sage/matrix/matrix_double_dense.pyx", line 1032:
    sage: sv = A.singular_values(eps='auto'); sv
Expected:
    verbose 1 (<module>) singular values, smallest-non-zero:cutoff:largest-zero, 2.2766...:6.2421...e-14:1.4160...e-15
    [35.139963659, 2.27661020871, 0.0, 0.0]
Got:
    verbose 1 (<module>) singular values, smallest-non-zero:cutoff:largest-zero, 2.27661020871:6.2421114782e-14:8.24999265856e-16
    [35.139963659, 2.27661020871, 0.0, 0.0]
**********************************************************************

So they are not too bad - an order of magnitude or so on (typically) very small values.

System is not too unusual. At Sage Days so do not have another system to test.

Very newly installed 64-bit Ubuntu 11.10 on Intel i7 M620.

$ gcc -v
Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/usr/lib/x86_64-linux-gnu/gcc/x86_64-linux-gnu/4.5.2/lto-wrapper
Target: x86_64-linux-gnu
Configured with: ../src/configure -v --with-pkgversion='Ubuntu/Linaro 4.5.2-8ubuntu4' --with-bugurl=file:///usr/share/doc/gcc-4.5/README.Bugs --enable-languages=c,c++,fortran,objc,obj-c++ --prefix=/usr --program-suffix=-4.5 --enable-shared --enable-multiarch --with-multiarch-defaults=x86_64-linux-gnu --enable-linker-build-id --with-system-zlib --libexecdir=/usr/lib/x86_64-linux-gnu --without-included-gettext --enable-threads=posix --with-gxx-include-dir=/usr/include/c++/4.5 --libdir=/usr/lib/x86_64-linux-gnu --enable-nls --with-sysroot=/ --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --enable-plugin --enable-gold --enable-ld=default --with-plugin-ld=ld.gold --enable-objc-gc --disable-werror --with-arch-32=i686 --with-tune=generic --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu
Thread model: posix
gcc version 4.5.2 (Ubuntu/Linaro 4.5.2-8ubuntu4) 

Rob

comment:21 Changed 10 years ago by jdemeyer

  • Reviewers changed from Martin Raum to Martin Raum, Jeroen Demeyer
  • Status changed from needs_review to positive_review

noise patch looks okay

comment:22 Changed 10 years ago by rbeezer

  • Keywords sd32 added

comment:23 Changed 10 years ago by leif

#10837, on which this ticket depends, needs work...

comment:24 Changed 10 years ago by jdemeyer

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