#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 )
Adds an exact routine, and a numerical routine, to determine if a matrix is Hermitian.
Apply:
Attachments (8)
Change History (38)
Changed 9 years ago by
comment:1 Changed 9 years ago by
- Status changed from new to needs_review
Changed 9 years ago by
comment:2 Changed 9 years ago by
- Description modified (diff)
comment:3 follow-up: ↓ 6 Changed 9 years ago by
typo: tranpose.
comment:4 follow-ups: ↓ 5 ↓ 7 Changed 9 years ago by
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
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
comment:7 in reply to: ↑ 4 ; follow-up: ↓ 8 Changed 9 years ago by
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
comment:8 in reply to: ↑ 7 ; follow-up: ↓ 9 Changed 9 years ago by
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: ↓ 10 Changed 9 years ago by
- Status changed from needs_review to needs_work
Replying to flawrence:
Would it be possible to copy the
matrix_double_dense.pyx
is_hermitian
intomatrix_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
andis_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: ↓ 11 Changed 9 years ago by
Replying to rbeezer:
Replying to flawrence:
Would it be possible to copy the
matrix_double_dense.pyx
is_hermitian
intomatrix_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
- 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: ↓ 13 Changed 9 years ago by
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
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
Alglib interface: #10880
comment:16 follow-up: ↓ 17 Changed 9 years ago by
- 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 9 years ago by
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 9 years ago by
comment:18 Changed 9 years ago by
Version 4 patch is a rebase to allow this to depend on #11027 for Schur decompositions.
Changed 9 years ago by
comment:19 Changed 9 years ago by
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 9 years ago by
comment:20 Changed 9 years ago by
- 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 9 years ago by
comment:21 Changed 9 years ago by
- 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 range
s so I just went with those. I think I am done. Really.
comment:22 follow-up: ↓ 23 Changed 9 years ago by
- 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 9 years ago by
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.
Apply: trac_10848-hermitian-matrices-v6.patch
comment:24 Changed 8 years ago by
- Dependencies set to #11027
comment:25 follow-up: ↓ 26 Changed 8 years ago by
- 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
comment:27 Changed 8 years ago by
- Milestone set to sage-4.7.2
comment:28 Changed 8 years ago by
- Description modified (diff)
comment:29 Changed 8 years ago by
- 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
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.
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.