Opened 11 years ago

# Fix and upgrade Gram-Schmidt — at Version 23

Reported by: Owned by: rbeezer jason, was major sage-4.8 linear algebra sd32 Rob Beezer Martin Raum 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

Note: See TracTickets for help on using tickets.