Opened 9 years ago

Closed 8 years ago

#10791 closed defect (fixed)

Fix and upgrade Gram-Schmidt

Reported by: rbeezer Owned by: jason, was
Priority: major Milestone: sage-4.8
Component: linear algebra Keywords: sd32
Cc: Merged in: sage-4.8.alpha5
Authors: Rob Beezer Reviewers: Martin Raum, John Palmieri
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by 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

Attachments (8)

trac_10791-fix-gram-schmidt.patch (26.7 KB) - added by rbeezer 9 years ago.
trac_10791-gram-schmidt-doctests.patch (2.1 KB) - added by rbeezer 9 years ago.
trac_10791-gram-schmidt-doctests-v2.patch (2.6 KB) - added by rbeezer 9 years ago.
trac-10791-gram-schmidt-doctests-v2.patch (2.6 KB) - added by mraum 9 years ago.
A copy of trac_10791-gram-schmidt-doctests-v2.patch
trac_10791-fix-gram-schmidt-v2.patch (26.8 KB) - added by rbeezer 8 years ago.
trac_10791-fix-gram-schmidt-v3.patch (26.9 KB) - added by rbeezer 8 years ago.
Apply only this patch
trac_10791-gram-schmidt-qqbar-doctests.patch (3.3 KB) - added by rbeezer 8 years ago.
trac_10791-gram-schmidt-negative-zero.patch (2.4 KB) - added by rbeezer 8 years ago.

Download all attachments as: .zip

Change History (42)

Changed 9 years ago by rbeezer

comment:1 Changed 9 years ago by rbeezer

  • Authors set to Rob Beezer
  • Description modified (diff)
  • Status changed from new to needs_review

comment:2 Changed 9 years ago by kini

  • Description modified (diff)

fix patchbot comment

comment:3 Changed 9 years ago by mraum

  • Reviewers set to Martin Raum
  • Status changed from needs_review to needs_work

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

Changed 9 years ago by rbeezer

comment:4 Changed 9 years ago by rbeezer

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

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

"doctests" patch implements changes to fix all three.

Rob

comment:5 Changed 9 years ago by mraum

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

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

comment:6 Changed 9 years ago by jdemeyer

  • Status changed from positive_review to needs_work

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.

Changed 9 years ago by rbeezer

comment:7 Changed 9 years ago by rbeezer

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

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.

comment:8 follow-up: Changed 9 years ago by mraum

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

comment:9 in reply to: ↑ 8 Changed 9 years ago by 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.

comment:10 Changed 9 years ago by rbeezer

Feeding the patchbot:

Depends on #10536, #10683, #10794

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

Changed 9 years ago by mraum

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

comment:11 Changed 9 years ago by 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

comment:12 follow-up: Changed 9 years ago by jdemeyer

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

comment:13 in reply to: ↑ 12 Changed 9 years ago by rbeezer

  • Status changed from needs_review to needs_info

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

Changed 8 years ago by rbeezer

comment:14 Changed 8 years ago by rbeezer

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

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.

Changed 8 years ago by rbeezer

Apply only this patch

comment:15 Changed 8 years ago by rbeezer

  • Description modified (diff)

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

comment:16 Changed 8 years ago by mraum

  • Status changed from needs_review to positive_review

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.

comment:17 Changed 8 years ago by rbeezer

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

Rob

comment:18 Changed 8 years ago by jdemeyer

  • Milestone changed from sage-4.7.1 to sage-4.7.2

comment:19 Changed 8 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:20 Changed 8 years ago by jdemeyer

  • Merged in sage-4.7.2.alpha2 deleted
  • Resolution fixed deleted
  • Status changed from closed to new

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]
**********************************************************************

comment:21 Changed 8 years ago by jdemeyer

  • Status changed from new to needs_review

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]
**********************************************************************

comment:22 Changed 8 years ago by jdemeyer

  • Status changed from needs_review to needs_work

comment:23 Changed 8 years ago by rbeezer

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

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

comment:24 Changed 8 years ago by rbeezer

  • Keywords sd32 added

comment:25 Changed 8 years ago by rbeezer

  • Status changed from needs_info to needs_review

comment:26 follow-up: Changed 8 years ago by 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).

comment:27 follow-up: Changed 8 years ago by jhpalmieri

  • Status changed from needs_review to needs_work

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

comment:28 in reply to: ↑ 26 Changed 8 years ago by 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.

comment:29 in reply to: ↑ 27 Changed 8 years ago by 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

Changed 8 years ago by rbeezer

comment:30 Changed 8 years ago by rbeezer

  • Description modified (diff)
  • Reviewers changed from Martin Raum to Martin Raum, John Palmieri
  • Status changed from needs_work to needs_review

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

comment:31 Changed 8 years ago by jason

  • Status changed from needs_review to positive_review

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)

comment:32 Changed 8 years ago by jdemeyer

  • Status changed from positive_review to needs_work

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

Changed 8 years ago by rbeezer

comment:33 Changed 8 years ago by rbeezer

  • Status changed from needs_work to positive_review

Jeroen,

My mistake - very sorry about that.

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

Rob

comment:34 Changed 8 years ago by jdemeyer

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