Opened 11 years ago

Last modified 10 years ago

#10791 closed defect

Fix and upgrade Gram-Schmidt — at Version 14

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

Status badges

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-v2.patch
  2. trac_10791-gram-schmidt-doctests-v2.patch

Depends on:

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

Change History (20)

Changed 11 years ago by rbeezer

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

Changed 11 years ago by rbeezer

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

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

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

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

Note: See TracTickets for help on using tickets.