Fix and upgrade Gram-Schmidt

[rbeezer]

If a set of vectors has a linearly dependent subset, followed by a vector outside the span of that subset, then the presence of a zero vector in the orthogonal set causes division by a zero dot product.

There are various other improvements that can be made in this routine.

sage: A=matrix(QQ, [[1,2],[2,4],[1,1]])
sage: A.gram_schmidt()
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)

/sage/sage-4.6.2.alpha4/<ipython console> in <module>()

/sage/sage-4.6.2.alpha4/local/lib/python2.6/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix.gram_schmidt (sage/matrix/matrix2.c:31979)()

/sage/sage-4.6.2.alpha4/local/lib/python2.6/site-packages/sage/modules/misc.pyc in gram_schmidt(B)
     55     for i in range(1,n):
     56         for j in range(i):
---> 57             mu[i,j] = B[i].dot_product(Bstar[j]) / (Bstar[j].dot_product(Bstar[j]))
     58         Bstar.append(B[i] - sum(mu[i,j]*Bstar[j] for j in range(i)))
     59     return Bstar, mu

/sage/sage-4.6.2.alpha4/local/lib/python2.6/site-packages/sage/structure/element.so in sage.structure.element.RingElement.__div__ (sage/structure/element.c:11981)()

/sage/sage-4.6.2.alpha4/local/lib/python2.6/site-packages/sage/rings/rational.so in sage.rings.rational.Rational._div_ (sage/rings/rational.c:15177)()

ZeroDivisionError: Rational division by zero

Patch summary:

• Original Gram-Schmidt has been adjusted to raise an error on linearly dependent input.

• Free module version of original Gram-Schmidt was used just one place, and it has been referenced directly, rather than through the matrix version.

• Matrix version is totally rewritten, building on QR factorizations as possible to get orthonormal bases. Linearly dependent subsets should be handled properly now. Return format has, by necessity, changed.

• The "noscale" helper method is reminiscent of the old code, but written in a column-oriented QR style.

Apply:

1. trac_10791-fix-gram-schmidt-v3.patch​
2. trac_10791-gram-schmidt-negative-zero.patch​
3. trac_10791-gram-schmidt-qqbar-doctests.patch​

Depends on:

1. #10536
2. #10683
3. #10794

[kini]

fix patchbot comment

[mraum]

In line 6364: properly contains In line 6444: For the integers and the rationals the field

There is only one thing you didn't catch. If you fix that I can give this a positive review (please provide a patch that depends on the current one, to make my life easier) :

File "/scratch/mraum/sagedev/devel/sage-dev/sage/modules/misc.py", line 66:

sage: gram_schmidt(V)

Exception raised:

Traceback (most recent call last):

File "/scratch/mraum/sagedev/local/bin/ncadoctest.py", line 1231, in run_one_test

self.run_one_example(test, example, filename, compileflags)

File "/scratch/mraum/sagedev/local/bin/sagedoctest.py", line 38, in run_one_example

OrigDocTestRunner?.run_one_example(self, test, example, filename, compileflags)

File "/scratch/mraum/sagedev/local/bin/ncadoctest.py", line 1172, in run_one_example

compileflags, 1) in test.globs

File "<doctest main.example_1[16]>", line 1, in <module>

gram_schmidt(V)###line 66:

sage: gram_schmidt(V)

File "/scratch/mraum/sagedev/local/lib/python/site-packages/sage/modules/misc.py", line 83, in gram_schmidt

raise ValueError?("linearly dependent input for module version of Gram-Schmidt")

ValueError?: linearly dependent input for module version of Gram-Schmidt

[rbeezer]

Martin - thanks for the careful review and for catching these three items.

"doctests" patch implements changes to fix all three.

Rob

[mraum]

Ok, the doctest patch fixes all issues. So this gets a positive review

[jdemeyer]

Problems while building the documentation:

docstring of sage.matrix.matrix2.Matrix.gram_schmidt:47: (WARNING/2) Literal block expected; none found.
docstring of sage.matrix.matrix2.Matrix.gram_schmidt:217: (WARNING/2) Literal block expected; none found.

[rbeezer]

Documentation fixes are added into a v2 version of the doctest patch.

Just had to change two :: to : in sage/matrix/matrix2.pyx.

So this change is ready for review.

[mraum]

The patchbot does not apply the patches this depends on. Do you, Robert, know how to fix that in the ticket description?

[rbeezer]

Replying to mraum:

The patchbot does not apply the patches this depends on. Do you, Robert, know how to fix that in the ticket description?

I have tried to chase back the dependencies and get everything in order for the patchbot. We'll see. The only change since this had a positive review was a couple of extra double-colons that had to be removed, so the documentation would build right.

[rbeezer]

Feeding the patchbot:

Depends on #10536, #10683, #10794

Apply: trac_10791-fix-gram-schmidt.patch, trac_10791-gram-schmidt-doctests-v2.patch

[mraum]

A copy of trac_10791-gram-schmidt-doctests-v2.patch

[mraum]

One more try to make the patchbot retry: I added a copy of the last patch.

Depends on #10536, #10683, #10794

Apply: trac_10791-fix-gram-schmidt.patch, trac-10791-gram-schmidt-doctests-v2.patch

[jdemeyer]

Using a bare except: (line 6303 of sage/matrix/matrix2.pyx) is a very bad idea. You should catch specific exceptions, not all.

[rbeezer]

Replying to jdemeyer:

Using a bare except: (line 6303 of sage/matrix/matrix2.pyx) is a very bad idea. You should catch specific exceptions, not all.

Yes, I became aware of this many months ago. However, I never know just what to be catching. In other words, it would be nice to use a __call___ method as a way to sanitize input, but I often have trouble figuring out just how the code will react when it fails (just what type of error will come back).

So which is preferable:

try:
   n = Integer(n_input)
except SomeError:
   raise ....

or

if not is_Integer(n_input):
    raise ...
else:
    n = Integer(n_input)

In the first case is it enough to just catch TypeError and ValueError?

In the second case, I am often at a loss to decide which to use, in particular I once spent a very long time trying to determine if an input was a real number.

So my naked except is really an act of desperation. Advice would be appreciated.

Rob

[rbeezer]

The naked except has been corrected to:

        try:
            F = R.fraction_field()
        except TypeError:
            raise TypeError(<snipped>)

This had a positive review back at comment 5. Since then it has received a documentation fix (double colons) and the fix on the except clause. Those are the only changes since then.

I do not know what is up with the patchbot, as it seems totally dead at the moment. In any event, it seems counter-productive to wait on it while it is still being developed. This is passing all tests for me right now.

[rbeezer]

Apply only this patch

[rbeezer]

Got confused with the various patches here, so I just rolled everything into one, the v3 patch. So only apply this one patch.

[mraum]

I checked this once more. The documentation builds fine, indeed. And everything since nothing but the exception has changed, I can set this again to positive review.

[rbeezer]

Thanks, Martin! I'm just getting caught up now that email notifications are working again.

Rob

[jdemeyer]

Doctest failure on flavius ():

sage -t -long  -force_lib devel/sage/sage/matrix/matrix2.pyx
**********************************************************************
File "/home/buildbot/build/sage/flavius-1/flavius_full/build/sage-4.7.2.alpha2/devel/sage-main/sage/matrix/matrix2.pyx", line 7722:
    sage: M.round(10)
Expected:
    [-2.4494897428           0.0           0.0]
    [-3.6742346142  0.7071067812           0.0]
    [-4.8989794856  1.4142135624           0.0]
Got:
    [-2.4494897428           0.0           0.0]
    [-3.6742346142  0.7071067812           0.0]
    [-4.8989794856  1.4142135624          -0.0]
**********************************************************************

[jdemeyer]

Doctest failure on bsd (OS X 10.8.0 x86_64):

sage -t -long  -force_lib devel/sage/sage/matrix/matrix2.pyx
**********************************************************************
File "/Users/buildbot/build/sage/bsd-1/bsd_64_full/build/sage-4.7.2.alpha2/devel/sage-main/sage/matrix/matrix2.pyx", line 7718:
    sage: G.round(6)
Expected:
    [-0.408248 -0.408248 -0.816497]
    [ 0.707107 -0.707107      -0.0]
    [ -0.57735  -0.57735   0.57735]
Got:
    [-0.408248 -0.408248 -0.816497]
    [ 0.707107 -0.707107       0.0]
    [ -0.57735  -0.57735   0.57735]
**********************************************************************

[rbeezer]

Turns out the IEEE specification includes a "negative zero." Rounding and .zero_at() might do the trick.

Six such changes in the "negative zero" patch. Successful with one instance where I converted a -0.0 to 0.0.

Jeroen - could you easily test/review this on the two failing machines? If it works across three architectures, I'll then apply the same technique to #10795, #10837 with some confidence.

Thanks for your patience with this (exasperating) numerical work.

Rob

[jason]

What's needed now? Jeroen just needs to test this last doctest fix on flavius and bsd? (or should we set it back to positive review, since the doctest changes are minor, and depend on the more widespread release testing to catch it on the platforms we don't have easy access to?)

I think I would prefer setting this back to positive review, since Jeroen will be implicitly testing it on those two platforms when he merges it (as long as it still passes on Rob's platform).

[jhpalmieri]

There are errors on sage.math and on OS X 10.6.8, applied to Sage 4.8.alpha3:

sage -t --long "devel/sage/sage/matrix/matrix2.pyx"         
**********************************************************************
File "/Applications/sage_builds/sage-4.8.alpha3/devel/sage/sage/matrix/matrix2.pyx", line 8072:
    sage: G*G.transpose()
Expected:
    [      1 0.?e-37 0.?e-35]
    [0.?e-37       1 0.?e-35]
    [0.?e-35 0.?e-35       1]
Got:
    [      1 0.?e-18 0.?e-19]
    [0.?e-18       1 0.?e-18]
    [0.?e-19 0.?e-18       1]
**********************************************************************
File "/Applications/sage_builds/sage-4.8.alpha3/devel/sage/sage/matrix/matrix2.pyx", line 8087:
    sage: G
Expected:
    [                         -0.3849001794597505?  -0.1924500897298753? - 0.1924500897298753?*I   0.3849001794597505? + 0.7698003589195010?*I                          -0.1924500897298753?]
    [-0.06165497274852388? - 0.1387236886841787?*I  0.8188551068163327? - 0.10018933071635130?*I  0.2003786614327026? + 0.05394810115495839?*I  0.02119389688230508? + 0.5028733714801478?*I]
    [  0.3842387256410419? - 0.5694103019142261?*I  0.1416892863096208? - 0.06139779741542298?*I  0.4633778333528464? - 0.01285039016180503?*I  0.02658516588219101? - 0.5373044261814995?*I]
Got:
    [                         -0.3849001794597505?  -0.1924500897298753? - 0.1924500897298753?*I   0.3849001794597505? + 0.7698003589195010?*I                          -0.1924500897298753?]
    [-0.06165497274852388? - 0.1387236886841787?*I  0.8188551068163327? - 0.10018933071635130?*I  0.2003786614327026? + 0.05394810115495839?*I  0.02119389688230508? + 0.5028733714801478?*I]
    [  0.3842387256410419? - 0.5694103019142261?*I   0.1416892863096208? - 0.0613977974154230?*I  0.4633778333528464? - 0.01285039016180503?*I  0.02658516588219101? - 0.5373044261814995?*I]
**********************************************************************
File "/Applications/sage_builds/sage-4.8.alpha3/devel/sage/sage/matrix/matrix2.pyx", line 8091:
    sage: M
Expected:
    [                        5.196152422706632?                                          0                                          0]
    [-3.079201435678004? + 3.464101615137755?*I             19.22286447225071? + 0.?e-37*I                                          0]
    [  4.426352063787131? - 8.66025403784439?*I -5.117362738127481? - 3.502773139275513?*I             7.480012456446966? + 0.?e-35*I]
Got:
    [                        5.196152422706632?                                          0                                          0]
    [-3.079201435678004? + 3.464101615137755?*I             19.22286447225071? + 0.?e-17*I                                          0]
    [  4.426352063787131? - 8.66025403784439?*I -5.117362738127481? - 3.502773139275513?*I             7.480012456446966? + 0.?e-16*I]
**********************************************************************
File "/Applications/sage_builds/sage-4.8.alpha3/devel/sage/sage/matrix/matrix2.pyx", line 8095:
    sage: M*G-A
Expected:
    [            0.?e-37 0.?e-37 + 0.?e-37*I 0.?e-37 + 0.?e-37*I             0.?e-37]
    [0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-36*I 0.?e-36 + 0.?e-36*I 0.?e-36 + 0.?e-36*I]
    [0.?e-35 + 0.?e-35*I 0.?e-35 + 0.?e-35*I 0.?e-35 + 0.?e-35*I 0.?e-35 + 0.?e-35*I]
Got:
    [            0.?e-18 0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-17*I             0.?e-18]
    [0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-17*I 0.?e-16 + 0.?e-17*I 0.?e-17 + 0.?e-16*I]
    [0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
**********************************************************************
File "/Applications/sage_builds/sage-4.8.alpha3/devel/sage/sage/matrix/matrix2.pyx", line 8099:
    sage: G*G.conjugate().transpose()
Expected:
    [1.000000000000000? + 0.?e-37*I            0.?e-37 + 0.?e-37*I            0.?e-36 + 0.?e-36*I]
    [           0.?e-37 + 0.?e-37*I 1.000000000000000? + 0.?e-37*I            0.?e-36 + 0.?e-36*I]
    [           0.?e-36 + 0.?e-36*I            0.?e-36 + 0.?e-36*I 1.000000000000000? + 0.?e-35*I]
Got:
    [1.000000000000000? + 0.?e-18*I            0.?e-18 + 0.?e-18*I            0.?e-17 + 0.?e-16*I]
    [           0.?e-18 + 0.?e-18*I 1.000000000000000? + 0.?e-18*I            0.?e-16 + 0.?e-17*I]
    [           0.?e-17 + 0.?e-16*I            0.?e-16 + 0.?e-17*I 1.000000000000000? + 0.?e-16*I]
**********************************************************************
File "/Applications/sage_builds/sage-4.8.alpha3/devel/sage/sage/matrix/matrix2.pyx", line 8130:
    sage: G*G.transpose()
Expected:
    [      1 0.?e-35]
    [0.?e-35       1]
Got:
    [      1 0.?e-18]
    [0.?e-18       1]
**********************************************************************
1 items had failures:
   6 of  68 in __main__.example_97
***Test Failed*** 6 failures.
For whitespace errors, see the file /Users/palmieri/.sage//tmp/matrix2_40831.py
	 [37.9 s]
 
----------------------------------------------------------------------
The following tests failed:


	sage -t --long "devel/sage/sage/matrix/matrix2.pyx"
Total time for all tests: 37.9 seconds

[jdemeyer]

Replying to jason:

and depend on the more widespread release testing to catch it on the platforms we don't have easy access to?

Yes. If you think a patch looks reasonable and has been tested to some extent, just set it to positive_review.

[rbeezer]

Replying to jhpalmieri:

There are errors on sage.math and on OS X 10.6.8, applied to Sage 4.8.alpha3:

Jason, John, Jeroen - thanks for the attention to this, and the other outstanding linear algebra tickets. I've been on an (extended) Sage vacation, but I think you are reeling me back in. ;-)

This does need work. Fix at #11544 prevents printing from having a side-effect that increases precision over QQbar. I think this work predates that change. So it does need work at the moment.

I am building a fresh alpha right now and will try to get on top of outstanding tickets later today. It'd be great to get a few of thse into 4.8 if it is not too late - I'll be using this stuff in class in the spring with my Sage-enabled textbook and a few of us are working on this CRC linear algebra handbook chapter that will be doctestable.

Thanks, Rob

[rbeezer]

Oh man, that was painful. I am really rusty. Use it or lose it. But I'm back.

Patch fixes failing doctests over QQbar, mostly by checking small norms on matrices which should be the zero matrix. A bit brutal, but we don't have round() and zero_at() for matrices over QQbar. (Those'd be a nice additions, no?).

I got identical failures on Ubuntu. Now passing all tests in matrix/matrix2.pyx. My 4.8.alpha3 has some weirdness (Atlas segfaults), but I think this should be OK.

Thanks for the look, John and Jason. Off to investigate a few other things.

Rob

[jason]

Looks good. Tested on OSX 10.6.8 and passes tests in matrix/*.pyx on sage-4.7.2.alpha3 (+some other patches I've applied)

[jdemeyer]

The commit message of trac_10791-gram-schmidt-negative-zero.patch​ needs to be updated.

[rbeezer]

Jeroen,

My mistake - very sorry about that.

I have replaced the problem patch with one containing a commit message.

Rob