sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.44k stars 481 forks source link

Checks for Hermitian matrices #10848

Closed 1d7ec08f-60ae-4512-91a6-8324c06eab9f closed 13 years ago

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago

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

Apply:

  1. attachment: trac_10848-hermitian-matrices-v7.patch

Depends on #11027

Component: linear algebra

Author: Rob Beezer

Reviewer: Mike Hansen

Merged: sage-4.7.2.alpha1

Issue created by migration from https://trac.sagemath.org/ticket/10848

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago

Author: Rob Beezer

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:1

Attachment: trac_10848-hermitian-matrices.patch.gz

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago

Description changed:

--- 
+++ 
@@ -1 +1,3 @@
 Adds an exact routine, and a numerical routine, to determine if a matrix is Hermitian.
+
+Apply trac_10848-hermitian-matrices-v2.patch
1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:2

Attachment: trac_10848-hermitian-matrices-v2.patch.gz

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.

56048686-9665-4c5e-973e-6c3add3aa805 commented 13 years ago
comment:3

typo: tranpose.

jasongrout commented 13 years ago
comment:4

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

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:5

Replying to @jasongrout:

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

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:6

Replying to @sagetrac-dsm:

typo: tranpose.

Thanks, got both of them fixed in latest patch.

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:7

Replying to @jasongrout:

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.
e2b0bcab-5fab-4c24-a972-fb341d16e224 commented 13 years ago
comment:8

Attachment: trac_10848-hermitian-matrices-v3.patch.gz

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?

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:9

Replying to @sagetrac-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.

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

jasongrout commented 13 years ago
comment:10

Replying to @rbeezer:

Replying to @sagetrac-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.

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:11

Replying to @jasongrout:

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.

jasongrout commented 13 years ago
comment:12

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.

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:13

Replying to @jasongrout:

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.

jasongrout commented 13 years ago
comment:14

Alglib interface: #10880

kini commented 13 years ago
comment:15

fix patchbot comment

kini commented 13 years ago

Description changed:

--- 
+++ 
@@ -1,3 +1,4 @@
 Adds an exact routine, and a numerical routine, to determine if a matrix is Hermitian.

 Apply trac_10848-hermitian-matrices-v2.patch
+Depends on #10536
18d65770-dc1e-4813-89a9-4828e4aae4a9 commented 13 years ago
comment:16

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

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:17

Replying to @sagetrac-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.

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:18

Attachment: trac_10848-hermitian-matrices-v4.patch.gz

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

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago

Attachment: trac_10848-hermitian-two-speed.patch.gz

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:19

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.

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago

Attachment: trac_10848-hermitian-matrices-v5.patch.gz

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:20

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.

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago

Description changed:

--- 
+++ 
@@ -1,4 +1,5 @@
 Adds an exact routine, and a numerical routine, to determine if a matrix is Hermitian.

-Apply trac_10848-hermitian-matrices-v2.patch
-Depends on #10536
+Apply trac_10848-hermitian-matrices-v5.patch
+
+Depends on #10536, #11027
1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago

Attachment: trac_10848-hermitian-matrices-v6.patch.gz

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago

Description changed:

--- 
+++ 
@@ -1,5 +1,5 @@
 Adds an exact routine, and a numerical routine, to determine if a matrix is Hermitian.

-Apply trac_10848-hermitian-matrices-v5.patch
+Apply trac_10848-hermitian-matrices-v6.patch

 Depends on #10536, #11027
1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:21

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.

18d65770-dc1e-4813-89a9-4828e4aae4a9 commented 13 years ago
comment:22

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.

18d65770-dc1e-4813-89a9-4828e4aae4a9 commented 13 years ago

Description changed:

--- 
+++ 
@@ -1,5 +1,8 @@
 Adds an exact routine, and a numerical routine, to determine if a matrix is Hermitian.

-Apply trac_10848-hermitian-matrices-v6.patch
+**Apply**:
+1. [attachment: trac_10848-hermitian-matrices-v6.patch](https://github.com/sagemath/sage-prod/files/10652213/trac_10848-hermitian-matrices-v6.patch.gz)

-Depends on #10536, #11027
+**Depends on**:
+1. #10536
+2. #11027
1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:23

Replying to @sagetrac-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

fchapoton commented 13 years ago

Dependencies: #11027

mwhansen commented 13 years ago

Reviewer: Mike Hansen

mwhansen commented 13 years ago
comment:25

Looks good to me.

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago
comment:26

Replying to @mwhansen:

Looks good to me.

Thanks for the review, Mike!

1d7ec08f-60ae-4512-91a6-8324c06eab9f commented 13 years ago

Attachment: trac_10848-hermitian-matrices-v7.patch.gz

Fixed commit message, rebased to 4.7.1.alpha4

jdemeyer commented 13 years ago

Description changed:

--- 
+++ 
@@ -1,8 +1,4 @@
 Adds an exact routine, and a numerical routine, to determine if a matrix is Hermitian.

 **Apply**:
-1. [attachment: trac_10848-hermitian-matrices-v6.patch](https://github.com/sagemath/sage-prod/files/10652213/trac_10848-hermitian-matrices-v6.patch.gz)
-
-**Depends on**:
-1. #10536
-2. #11027
+1. [attachment: trac_10848-hermitian-matrices-v7.patch](https://github.com/sagemath/sage-prod/files/10652214/trac_10848-hermitian-matrices-v7.patch.gz)
jdemeyer commented 13 years ago

Merged: sage-4.7.2.alpha1

jdemeyer commented 7 years ago
comment:30

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.