Checks for Hermitian matrices

[rbeezer]

Adds an exact routine, and a numerical routine, to determine if a matrix is Hermitian.

Apply:

1. trac_10848-hermitian-matrices-v7.patch​

[rbeezer]

Had an off-by-one error and was not checking the diagonal elements. Fixed now in the v2 patch, and added a doctest that would have caught the mistake.

[dsm]

typo: tranpose.

[jason]

Can you clarify this statement? "For numerical matrices a specialized routine available over RDF and CDF is a good choice. " When I read it, I'm not sure what it means---should I program my own routine?

Maybe you could change it to: "For numerical matrices over RDF or CDF, the tolerance for comparison can also be specified (see ~REFERENCE)."

[rbeezer]

Replying to jason:

Can you clarify this statement? "For numerical matrices a specialized routine available over RDF and CDF is a good choice. " When I read it, I'm not sure what it means---should I program my own routine?

Yes, that could be improved (and I used the same thing on some other ticket). What I was trying to convey was the idea that RDF/CDF are better for numerical work than RR/CC. There are many methods designed for exact rings that get applied to RDF/CDF/RR/CC, and I am hoping to improve the situation for RDF/CDF, thus an effort to steer folks there. Maybe I should just say that outright, plus mention the tolerance option as you have suggested.

I'll have a new patch up later today.

Maybe you could change it to: "For numerical matrices over RDF or CDF, the tolerance for comparison can also be specified (see ~REFERENCE)."

[rbeezer]

Replying to dsm:

typo: tranpose.

Thanks, got both of them fixed in latest patch.

[rbeezer]

Replying to jason:

Can you clarify this statement?

Does this sound better? Let me know and I'll replicate into is_unitary.

        This routine is for matrices over exact rings and so may not
        work properly for matrices over ``RR`` or ``CC``.  For matrices with
        approximate entries, the rings of double-precision floating-point
        numbers, ``RDF`` and ``CDF``, are a better choice since the
        :meth:`sage.matrix.matrix_double_dense.Matrix_double_dense.is_hermitian`
        method has a tolerance parameter.  This provides control over
        allowing for minor discrepancies between entries when checking
        equality.

[flawrence]

Replying to rbeezer:

        This routine is for matrices over exact rings and so may not
        work properly for matrices over ``RR`` or ``CC``.  For matrices with
        approximate entries, the rings of double-precision floating-point
        numbers, ``RDF`` and ``CDF``, are a better choice since the
        :meth:`sage.matrix.matrix_double_dense.Matrix_double_dense.is_hermitian`
        method has a tolerance parameter.  This provides control over
        allowing for minor discrepancies between entries when checking
        equality.

Would it be possible to copy the matrix_double_dense.pyx is_hermitian into matrix_dense.pyx (adjusting the default tolerance and doctests, of course) and thus remove the quirk that makes this warning necessary?

[rbeezer]

Replying to flawrence:

Would it be possible to copy the matrix_double_dense.pyx is_hermitian into matrix_dense.pyx (adjusting the default tolerance and doctests, of course) and thus remove the quirk that makes this warning necessary?

Hi Felix,

This sounds like a good idea.

• I'd imagine the code in matrix2 branching for exact vs. inexact rings. Any tolerance would be ignored for exact rings.

• Are RDF/CDF/RR/CC the only inexact rings in Sage? They need to be amenable to an absolute value in order to do the comparison. As organized in the patch, we at least know just which ring we are dealing with.

• Same idea would apply to is_symmetric and is_unitary.

• I'm guessing this will change some behavior if applied to is_symmetric. In other words, I bet the exact version gets called for some inexact rings. I may test this later.

BTW, I saw your post on sage-devel asking for greater SciPY/NumPy integration. I'm hoping to (slowly) make more of the matrix algebra available, so maybe that will help.

Jason - any comments on the above?

Rob

[jason]

Replying to rbeezer:

Replying to flawrence:

Would it be possible to copy the matrix_double_dense.pyx is_hermitian into matrix_dense.pyx (adjusting the default tolerance and doctests, of course) and thus remove the quirk that makes this warning necessary?

Hi Felix,

This sounds like a good idea.

• I'd imagine the code in matrix2 branching for exact vs. inexact rings. Any tolerance would be ignored for exact rings.

If the parameter is there, don't ignore it. That would be really confusing.

• Are RDF/CDF/RR/CC the only inexact rings in Sage? They need to be amenable to an absolute value in order to do the comparison. As organized in the patch, we at least know just which ring we are dealing with.

sage: SR.is_exact()
False

In fact, there are an infinite number of inexact rings:

sage: RealField(100).is_exact()
False
sage: S.<s> = LaurentSeriesRing(GF(5))
sage: T.<t> = PowerSeriesRing(pAdicRing(5))
sage: T.is_exact()
False

There probably many, many more.

[rbeezer]

Replying to jason:

If the parameter is there, don't ignore it. That would be really confusing.

I wouldn't totally ignore it, but I'd not honor it either, I think. Throw an error, since it would exhibit a basic misunderstanding/misapplication.

sage: SR.is_exact()
False

That's the one that would bite me. I knew there was one oddball one.

In fact, there are an infinite number of inexact rings:

sage: RealField(100).is_exact()
False
sage: S.<s> = LaurentSeriesRing(GF(5))
sage: T.<t> = PowerSeriesRing(pAdicRing(5))
sage: T.is_exact()
False

Thanks, that's what I needed to know. (I use RR/CC as stand-ins for all the RealField()'s.)

So this will be a problem:

sage: T.<t> = PowerSeriesRing(pAdicRing(5))
sage: a=T.an_element()
sage: a
(1 + O(5^20))*t
sage: abs(a)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

/sage/dev/devel/sage-main/<ipython console> in <module>()

TypeError: bad operand type for abs(): 'sage.rings.power_series_poly.PowerSeries_poly'

Which is making me think it would be better to make the annoying message go away by "somebody" implementing matrices over RR/CC/RealField()/ComplexField?() properly.

[jason]

Yep, I'm really looking forward to an alglib interface, which seems like the best contender right now for a good RealField/ComplexField? matrix class backend.

[rbeezer]

Replying to jason:

Yep, I'm really looking forward to an alglib interface, which seems like the best contender right now for a good RealField/ComplexField? matrix class backend.

Aah, that looks nice.

[jason]

Alglib interface: #10880

[kini]

fix patchbot comment

[mraum]

I suggest to reword line 2945 like "A matrix that is nearly Hermitian, but for one non-real" and I would introduce one keyword for the RDF implementation deciding whether the entries are naively compared (quick and what students might assume) or the svd is applied and the imaginary values considered (numerically much better conditioned).

[rbeezer]

Replying to mraum:

I suggest to reword line 2945 like "A matrix that is nearly Hermitian, but for one non-real"

Thanks, Martin. I'll make that change.

So a check based on the Schur decomposition at #11027 will be a good high-reliability test, while the naive cut-off comparison can be a high-speed crude check.

[rbeezer]

Version 4 patch is a rebase to allow this to depend on #11027 for Schur decompositions.

[rbeezer]

Now depends on #11027, then apply v4 patch, then "two-speed" patch.

This is not ready, mostly posted for safe-keeping. Needs more docs, maybe some timing tests. But it should work and only one test fails (needs tolerance adjustment, I'd guess). More soon.

[rbeezer]

v5 patch is self-contained, apply only this one.

Two options for the check, the naive one, or one based on the Schur decomposition (#11027).

This needs to check that the upper half of a matrix is zero, so I broke out a helper method for that, since I'll use it in a future is_normal() method. I might pair it with an upper-triangular check at some point and make them both visible. But not as part of this.

[rbeezer]

v6 adds caching the exact version's result, and fixes some (more) off-by-one problems with the old-style loops. Cython now efficiently translates ranges so I just went with those. I think I am done. Really.

[mraum]

For some reason the patchbot does not apply this correctly. The changes to the description should fix this. If so I will come back to this.

[rbeezer]

Replying to mraum:

For some reason the patchbot does not apply this correctly. The changes to the description should fix this. If so I will come back to this.

I think the description may be for the release manager, and the buildbot reads comments.

​http://wiki.sagemath.org/buildbot

Also, the buildbot is "stuck" back at 4.6.2. Maybe the following will help.

Depends on #10536, #11027

Apply: trac_10848-hermitian-matrices-v6.patch

[mhansen]

Looks good to me.

[rbeezer]

Replying to mhansen:

Looks good to me.

Thanks for the review, Mike!

[rbeezer]

Fixed commit message, rebased to 4.7.1.alpha4

[jdemeyer]

Does anybody here remember why numpy.absolute(numpy.imag(z)) is called on a CDF element z instead of simply abs(z.imag())? The former breaks with numpy 1.13 and I'm changing the code to the latter on #24063.