Opened 9 years ago

Closed 8 years ago

Last modified 2 years ago

#10848 closed enhancement (fixed)

Checks for Hermitian matrices

Reported by: rbeezer Owned by: jason, was
Priority: minor Milestone: sage-4.7.2
Component: linear algebra Keywords:
Cc: Merged in: sage-4.7.2.alpha1
Authors: Rob Beezer Reviewers: Mike Hansen
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: #11027 Stopgaps:

Description (last modified by jdemeyer)

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

Apply:

  1. trac_10848-hermitian-matrices-v7.patch

Attachments (8)

trac_10848-hermitian-matrices.patch (6.0 KB) - added by rbeezer 9 years ago.
trac_10848-hermitian-matrices-v2.patch (6.8 KB) - added by rbeezer 9 years ago.
trac_10848-hermitian-matrices-v3.patch (7.3 KB) - added by rbeezer 9 years ago.
trac_10848-hermitian-matrices-v4.patch (7.0 KB) - added by rbeezer 8 years ago.
trac_10848-hermitian-two-speed.patch (5.0 KB) - added by rbeezer 8 years ago.
trac_10848-hermitian-matrices-v5.patch (11.0 KB) - added by rbeezer 8 years ago.
trac_10848-hermitian-matrices-v6.patch (11.8 KB) - added by rbeezer 8 years ago.
trac_10848-hermitian-matrices-v7.patch (11.6 KB) - added by rbeezer 8 years ago.
Fixed commit message, rebased to 4.7.1.alpha4

Download all attachments as: .zip

Change History (38)

Changed 9 years ago by rbeezer

comment:1 Changed 9 years ago by rbeezer

  • Authors set to Rob Beezer
  • Status changed from new to needs_review

Changed 9 years ago by rbeezer

comment:2 Changed 9 years ago by rbeezer

  • Description modified (diff)

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.

comment:3 follow-up: Changed 9 years ago by dsm

typo: tranpose.

comment:4 follow-ups: Changed 9 years ago by 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)."

comment:5 in reply to: ↑ 4 Changed 9 years ago by 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)."

comment:6 in reply to: ↑ 3 Changed 9 years ago by rbeezer

Replying to dsm:

typo: tranpose.

Thanks, got both of them fixed in latest patch.

comment:7 in reply to: ↑ 4 ; follow-up: Changed 9 years ago by 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.

Changed 9 years ago by rbeezer

comment:8 in reply to: ↑ 7 ; follow-up: Changed 9 years ago by 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?

comment:9 in reply to: ↑ 8 ; follow-up: Changed 9 years ago by rbeezer

  • Status changed from needs_review to needs_work

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

comment:10 in reply to: ↑ 9 ; follow-up: Changed 9 years ago by 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.

comment:11 in reply to: ↑ 10 Changed 9 years ago by rbeezer

  • Status changed from needs_work to needs_review

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.

comment:12 follow-up: Changed 9 years ago by 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.

comment:13 in reply to: ↑ 12 Changed 9 years ago by 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.

comment:14 Changed 9 years ago by jason

Alglib interface: #10880

comment:15 Changed 8 years ago by kini

  • Description modified (diff)

fix patchbot comment

comment:16 follow-up: Changed 8 years ago by mraum

  • Status changed from needs_review to needs_work

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).

comment:17 in reply to: ↑ 16 Changed 8 years ago by 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.

Changed 8 years ago by rbeezer

comment:18 Changed 8 years ago by rbeezer

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

Changed 8 years ago by rbeezer

comment:19 Changed 8 years ago by 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.

Changed 8 years ago by rbeezer

comment:20 Changed 8 years ago by rbeezer

  • Description modified (diff)
  • Status changed from needs_work to needs_review

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.

Changed 8 years ago by rbeezer

comment:21 Changed 8 years ago by rbeezer

  • Description modified (diff)

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.

comment:22 follow-up: Changed 8 years ago by mraum

  • Description modified (diff)

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.

comment:23 in reply to: ↑ 22 Changed 8 years ago by 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

comment:24 Changed 8 years ago by chapoton

  • Dependencies set to #11027

comment:25 follow-up: Changed 8 years ago by mhansen

  • Reviewers set to Mike Hansen
  • Status changed from needs_review to positive_review

Looks good to me.

comment:26 in reply to: ↑ 25 Changed 8 years ago by rbeezer

Replying to mhansen:

Looks good to me.

Thanks for the review, Mike!

comment:27 Changed 8 years ago by jdemeyer

  • Milestone set to sage-4.7.2

Changed 8 years ago by rbeezer

Fixed commit message, rebased to 4.7.1.alpha4

comment:28 Changed 8 years ago by jdemeyer

  • Description modified (diff)

comment:29 Changed 8 years ago by jdemeyer

  • Merged in set to sage-4.7.2.alpha1
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:30 Changed 2 years ago by 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.

Note: See TracTickets for help on using tickets.