Opened 10 years ago
Closed 9 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 )
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:
- trac_10791-fix-gram-schmidt-v3.patch
- trac_10791-gram-schmidt-negative-zero.patch
- trac_10791-gram-schmidt-qqbar-doctests.patch
Depends on:
Attachments (8)
Change History (42)
Changed 10 years ago by
comment:1 Changed 10 years ago by
- Description modified (diff)
- Status changed from new to needs_review
comment:2 Changed 10 years ago by
- Description modified (diff)
comment:3 Changed 10 years ago by
- 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 10 years ago by
comment:4 Changed 10 years ago by
- 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
- 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
- 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 10 years ago by
comment:7 Changed 10 years ago by
- 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
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
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 10 years ago by
comment:11 Changed 10 years ago by
comment:12 follow-up: ↓ 13 Changed 10 years ago by
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
- Status changed from needs_review to needs_info
Replying to jdemeyer:
Using a bare
except:
(line 6303 ofsage/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
comment:14 Changed 10 years ago by
- 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.
comment:15 Changed 10 years ago by
- 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
- 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
Thanks, Martin! I'm just getting caught up now that email notifications are working again.
Rob
comment:18 Changed 10 years ago by
- Milestone changed from sage-4.7.1 to sage-4.7.2
comment:19 Changed 10 years ago by
- 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
- 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
- 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
- Status changed from needs_review to needs_work
comment:23 Changed 10 years ago by
- 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 10 years ago by
- Keywords sd32 added
comment:25 Changed 10 years ago by
- Status changed from needs_info to needs_review
comment:26 follow-up: ↓ 28 Changed 9 years ago by
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 9 years ago by
- 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 9 years ago by
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 9 years ago by
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 9 years ago by
comment:30 Changed 9 years ago by
- 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 9 years ago by
- 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 9 years ago by
- Status changed from positive_review to needs_work
The commit message of trac_10791-gram-schmidt-negative-zero.patch needs to be updated.
Changed 9 years ago by
comment:33 Changed 9 years ago by
- 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 9 years ago by
- Merged in set to sage-4.8.alpha5
- Resolution set to fixed
- Status changed from positive_review to closed
fix patchbot comment