Singular values of matrices over CDF

[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

[mraum]

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

[jdemeyer]

This patch conflicts with #10837.

[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.}

[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.

[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

[rbeezer]

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.

[mraum]

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.

[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

[rbeezer]

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.

[jdemeyer]

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.

[rbeezer]

Rebase

[rbeezer]

Numerical noise in doctests

[rbeezer]

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.

[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.

[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

[jdemeyer]

noise patch looks okay

[leif]

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