#10791: Fix and upgrade Gram-Schmidt
-------------------------------+--------------------------------------------
Reporter: rbeezer | Owner: jason, was
Type: defect | Status: needs_info
Priority: major | Milestone: sage-4.7.2
Component: linear algebra | Resolution:
Keywords: | Work_issues:
Upstream: N/A | Reviewer: Martin Raum
Author: Rob Beezer | Merged:
Dependencies: |
-------------------------------+--------------------------------------------
Changes (by rbeezer):
* status: needs_work => needs_info
Old description:
> 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. [attachment:trac_10791-fix-gram-schmidt-v3.patch]
>
> '''Depends on:'''
> 1. #10536
> 2. #10683
> 3. #10794
New description:
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. [attachment:trac_10791-fix-gram-schmidt-v3.patch]
1. [attachment:trac_10791-gram-schmidt-negative-zero.patch]
'''Depends on:'''
1. #10536
2. #10683
3. #10794
--
Comment:
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
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/10791#comment:23>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
