Gram-Schmidt Procedure for Symbolic Ring

[Michael Jung]

Hello everyone,

I have a question regarding the implemented Gram-Schmidt procedure. Using a matrix in the symbolic ring leads to the following error message:

m = matrix(SR, [[x,2*x],[x^2,1]])
m.gram_schmidt()

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-21-6b8bbf989823> in <module>()
      1 m = matrix(SR, [[x,Integer(2)*x],[x**Integer(2),Integer(1)]])
----> 2 m.gram_schmidt()

/home/michi/GitProjects/sage/local/lib/python3.7/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.gram_schmidt (build/cythonized/sage/matrix/matrix2.c:70206)()
   9871                 Q, R = self.transpose()._gram_schmidt_noscale()
   9872         else:
-> 9873             raise NotImplementedError("Gram-Schmidt orthogonalization not implemented for matrices over inexact rings, except for RDF and CDF")
   9874         return Q.transpose(), R.transpose()
   9875 

NotImplementedError: Gram-Schmidt orthogonalization not implemented for matrices over inexact rings, except for RDF and CDF

Why isn't it possible to compute an orthonormal basis from vectors with entries in the symbolic ring? I see no problem in that. What happens behind the scenes?

Thanks in advance and best wishes

Michael

[Emmanuel Charpentier]

Well... the error message is pretty explicit: since

sage: SR.is_exact()

False

M.gram_schmidt() wont work if M.base_ring is SR.

Creating a second special case for SR may not be as simple as for RDF, since a lot of other cases (beyond RDF) can happen in this case....

[Michael Jung]

Maybe, I did get something wrong. But what's the problem about Gram-Schmidt on SR? There are just sums and divisions (and probably roots to normalize) in Gram-Schmidt which should not lead to problems in SR.

By the way, what does "exact" actually mean?

[Simon King]

Hi Michael,

On 2019-10-24, Michael Jung <mic...@uni-potsdam.de> wrote:
> Maybe, I did get something wrong. But what's the problem about Gram-Schmidt
> on SR? There are just sums and divisions (and probably roots to normalize)
> in Gram-Schmidt which should not lead to problems in SR.
>
> By the way, what does "exact" actually mean?

It means that arithmetics is not affected by rounding/approximations.
And some methods (apparently Gram-Schmidt is one of them) is only "safe"
when no rounding occurs.

That said, I'd appreciate to add an option that forces such methods
even in an inexact ring.

Best regards,
Simon

[Michael Jung]

Do you have an example where SR fails to be exact?

[Simon King]

On 2019-10-24, Michael Jung <mic...@uni-potsdam.de> wrote:
> Do you have an example where SR fails to be exact?

One can convert a float to SR. The result is in SR, but still behaves
like a float:
sage: a = SR(2.)^(1/500)
sage: a^500
2.00000000000005
sage: a.parent()
Symbolic Ring

Best regards,
Simon

[Michael Jung]

I see. Maybe it is possible to decompose/split the matrix in SR into an exact and inexact part, convert the inexact part to RDF and apply the Gram-Schmidt algorithm appropiately? I don't know, maybe it's too naive?

[Emmanuel Charpentier]

Dear Michael

Contrast

sage: A=matrix(QQ,[[1, 2], [3,4]])

sage: G,M=M.gram_schmidt()

sage: A=matrix(QQ,[[1, 2], [3,4]])

sage: G,M=A.gram_schmidt()

sage: M*G==A

True

with

sage: Ar=matrix(RDF,[[1, 2], [3,4]])

sage: Gr,Mr=Ar.gram_schmidt()

sage: Mr*Gr==Ar

False

sage: Mr*Gr-A

[-1.1102230246251565e-16  -2.220446049250313e-16]

[-1.3322676295501878e-15                     0.0]

Writing *correctly* this decomposition is, IIRC, a numerical analysis bitch... You are, IIRC, led to compute differences of large products, where underflows can easily slip into... Definitely not an amateur's problem.

While I agree in principle with Simon, actually *doing* it isn't easy.

[Simon King]

Hi Emmanuel,

On 2019-10-24, Emmanuel Charpentier <emanuel.c...@gmail.com> wrote:
> Writing *correctly* this decomposition is, IIRC, a numerical analysis
> bitch... You are, IIRC, led to compute differences of large products, where
> underflows can easily slip into... Definitely not an amateur's problem.
>
> While I agree in principle with Simon, actually *doing* it isn't easy.

Indeed. I had teaching in mind when I suggested to have an option to
force "numerically evil" methods, as the students tend to not believe
that a method they learned in the first year can be problematic.

Best regards,
Simon

[Vincent Delecroix]

This was an easy one. The following shows that SR is just
broken.... pi is rational!

sage: q = continued_fraction(pi).convergent(100)
sage: q
8736149038303113005348154524599771853409352442745266/2780802606066896232581239559281727773240004199722661
sage: bool(pi == q)
True

And everything is with exact numbers.

Vincent

[Simon King]

Hi Vincent,

On 2019-10-27, Vincent Delecroix <20100.d...@gmail.com> wrote:
> This was an easy one. The following shows that SR is just
> broken.... pi is rational!
>
> sage: q = continued_fraction(pi).convergent(100)
> sage: q
> 8736149038303113005348154524599771853409352442745266/2780802606066896232581239559281727773240004199722661
> sage: bool(pi == q)
> True

Nice!

Best regards,
Simon

[Michael Jung]

Nice example! I see your point.

However, I wonder. The matrix inversion should cause the same problem, but it is implemented for the symbolic ring. What is different?

Best wishes
Michael

[Nils Bruin]

On Monday, October 28, 2019 at 3:30:03 AM UTC-7, Michael Jung wrote:

Nice example! I see your point.

However, I wonder. The matrix inversion should cause the same problem, but it is implemented for the symbolic ring. What is different?

I don't think it is really implemented *for* the symbolic ring. It's just an implementation that gets used for the symbolic ring and once people submit enough bug reports about it, it probably needs to be deprecated of be adorned with a warning message.

For Gram-Schmidt, there's a further problem: you need a definite inner product for it to work. SR contains elements from C. So do we take the standard Hermitian inner product? In fact, there are all kinds of functions in there as well, so perhaps a Hilbert space inner product is appropriate. But which one? I think the problems of SR for this are so obvious that no-one dared to ignore them when thinking of providing an implementation.

[Michael Orlitzky]

On 10/28/19 12:03 PM, Nils Bruin wrote:
> On Monday, October 28, 2019 at 3:30:03 AM UTC-7, Michael Jung wrote:
>
> Nice example! I see your point.
>
> However, I wonder. The matrix inversion should cause the same
> problem, but it is implemented for the symbolic ring. What is different?
>
>
> I don't think it is really implemented *for* the symbolic ring. It's
> just an implementation that gets used for the symbolic ring and once
> people submit enough bug reports about it, it probably needs to be
> deprecated of be adorned with a warning message.
>
> For Gram-Schmidt, there's a further problem: you need a definite inner

> product for it to work....

I'm torn on this one. I think the practical reason that Gram-Schmidt
isn't implemented for inexact rings is because the QR method isn't
implemented for matrices over inexact rings, and all gram_schmidt() does
is call QR().

And while I know that linear algebra doesn't work over an inexact field,
it would still be nice sometimes if I didn't have to play these games in
my library code. Ultimately, all you need for Gram-Schmidt is a list of
vectors in some space where an inner_product() method is defined. No
matrix or anything, just a list of vectors. Then you compute the thing
using the Wikipedia algorithm, and get back another list of orthogonal
(and optionally normal) vectors.

There's no need to even involve a matrix, and having to go from an
abstract Hilbert space element to a vector to a matrix to the matrix's
transpose to its QR factorization to its transpose to a list of vectors
that might be zero to a list of vectors that are non-zero and then back
to abstract elements starts to feel like maybe it's not the most
efficient way to do things.

Will it give you great answers on pathological sets? Nope. But nothing
else in linear algebra works either, so as long as we're still allowed
to create vector spaces over RR, SR, and friends -- this seems like an
odd place to draw a line in the sand. I've had to reimplement
Gram-Schmidt four or five times now to work around this so that people
who use my code can ignore my warnings and use RealField(100) as their
base field.