Opened 10 years ago

Last modified 4 years ago

#10848 closed enhancement

Checks for Hermitian matrices — at Version 21

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

Status badges

Description (last modified by rbeezer)

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

Apply trac_10848-hermitian-matrices-v6.patch

Depends on #10536, #11027

Change History (28)

Changed 10 years ago by rbeezer

comment:1 Changed 10 years ago by rbeezer

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

Changed 10 years ago by rbeezer

comment:2 Changed 10 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 10 years ago by dsm

typo: tranpose.

comment:4 follow-ups: Changed 10 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 10 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 10 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 10 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 10 years ago by rbeezer

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

Alglib interface: #10880

comment:15 Changed 10 years ago by kini

  • Description modified (diff)

fix patchbot comment

comment:16 follow-up: Changed 10 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 10 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 10 years ago by rbeezer

comment:18 Changed 10 years ago by rbeezer

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

Changed 10 years ago by rbeezer

comment:19 Changed 10 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 10 years ago by rbeezer

comment:20 Changed 10 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 10 years ago by rbeezer

comment:21 Changed 10 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.

Note: See TracTickets for help on using tickets.