Opened 11 years ago

# Fix and upgrade Gram-Schmidt — at Version 30

Reported by: Owned by: rbeezer jason, was major sage-4.8 linear algebra sd32 Rob Beezer Martin Raum, John Palmieri N/A

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:

Depends on:

### comment:1 Changed 11 years ago by rbeezer

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

### comment:2 Changed 11 years ago by kini

• Description modified (diff)

fix patchbot comment

### comment:3 Changed 11 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

### comment:4 Changed 11 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 10 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 10 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.
```

### comment:7 Changed 10 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: ↓ 9 Changed 10 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 10 years ago by rbeezer

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 10 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 10 years ago by mraum

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

### comment:11 Changed 10 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: ↓ 13 Changed 10 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 10 years ago by rbeezer

• Status changed from needs_review to needs_info

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

### comment:14 Changed 10 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 10 years ago by rbeezer

Apply only this patch

### comment:15 Changed 10 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 10 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 10 years ago by rbeezer

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

Rob

### comment:18 Changed 10 years ago by jdemeyer

• Milestone changed from sage-4.7.1 to sage-4.7.2

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

• Status changed from needs_review to needs_work

### comment:23 Changed 10 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:25 Changed 10 years ago by rbeezer

• Status changed from needs_info to needs_review

### comment:26 follow-up: ↓ 28 Changed 10 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: ↓ 29 Changed 10 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]
**********************************************************************
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 10 years ago by jdemeyer

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 10 years ago by rbeezer

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

### comment:30 Changed 10 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

Note: See TracTickets for help on using tickets.