Opened 10 years ago
Last modified 10 years ago
#11027 closed enhancement
Schur matrix decomposition over RDF/CDF — at Version 20
Reported by: | rbeezer | Owned by: | jason, was |
---|---|---|---|
Priority: | minor | Milestone: | sage-4.7.1 |
Component: | linear algebra | Keywords: | |
Cc: | jason, kcrisman | Merged in: | |
Authors: | Rob Beezer | Reviewers: | Martin Raum |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
Wraps the schur()
function from SciPy
.
Apply:
Change History (22)
Changed 10 years ago by
comment:1 follow-up: ↓ 2 Changed 10 years ago by
- Reviewers set to Martin Raum
- Status changed from new to needs_work
comment:2 in reply to: ↑ 1 Changed 10 years ago by
Replying to mraum:
Some discussion is necessary
Definitely! That is disturbing. I don't know if SciPy
has a random element, or it is that platform dependent.
This decomposition is not unique (mathematically). I cannot determine if SciPy
is suppose to return the same result across platforms. I'll dig deeper.
In any event, your three examples all seem correct, in that the diagonals of T all have the eigenvalues of the matrix. Run in Sage, the doctests should have checked all the key properties (unitary, upper-triangular, decomposition).
Did you check Mathematica's result for unitary and decomposition?
Rob
comment:3 follow-up: ↓ 6 Changed 10 years ago by
A = matrix(RDF, [[-7., 5., 11., -4., 13.], [-11., -3., 11., 8., -19.], [-6., 3., -5., 0., -12.], [-4., -12., -14., 8., -8.], [11., 0., 9., 6., 10.]]) Q, T = A.schur(base_ring=CDF) T.round(4)
Unknown hardware, etc
http://abel.ee.ucla.edu/cvxopt/userguide/lapack.html
[ 5.67e+00+j1.69e+01 -2.13e+01+j2.85e+00 1.40e+00+j5.88e+00 -4.19e+00+j2.05e-01 3.19e+00-j1.01e+01] [ 0.00e+00-j0.00e+00 5.67e+00-j1.69e+01 1.09e+01+j5.93e-01 -3.29e+00-j1.26e+00 -1.26e+01+j7.80e+00] [ 0.00e+00-j0.00e+00 0.00e+00-j0.00e+00 1.27e+01+j3.43e-17 -6.83e+00+j2.18e+00 5.31e+00-j1.69e+00] [ 0.00e+00-j0.00e+00 0.00e+00-j0.00e+00 0.00e+00-j0.00e+00 -1.31e+01-j0.00e+00 -2.60e-01-j0.00e+00] [ 0.00e+00-j0.00e+00 0.00e+00-j0.00e+00 0.00e+00-j0.00e+00 0.00e+00-j0.00e+00 -7.86e+00-j0.00e+00]
sage-4.7.alpha2, sage.math
[ 5.6668 - 16.9373*I -8.6458 + 7.8734*I -7.1561 + 14.5703*I 5.7811 + 2.688*I -8.5999 - 5.1071*I] [ 0 5.6668 + 16.9373*I -10.0921 + 9.9632*I -1.5618 - 5.8776*I 7.4924 + 8.7084*I] [ 0 0 12.6587 1.6512 + 1.049*I 11.1778 + 4.4387*I] [ 0 0 0 -13.136 -0.1856 + 0.0353*I] [ 0 0 0 0 -7.8563]
sage4.7.alpha2, Intel i7-2600
[ 5.6668 - 16.9373*I -8.6458 + 7.8734*I -7.1561 + 14.5703*I 5.7811 + 2.688*I -8.5999 - 5.1071*I] [ 0 5.6668 + 16.9373*I -10.0921 + 9.9632*I -1.5618 - 5.8776*I 7.4924 + 8.7084*I] [ 0 0 12.6587 1.6512 + 1.049*I 11.1778 + 4.4387*I] [ 0 0 0 -13.136 -0.1856 + 0.0353*I] [ 0 0 0 0 -7.8563]
Mathematica 6.0, sage.math
{ {5.66679 - 16.9373 I, -10.6919 + 4.73528 I, -7.15609 + 14.5703 I, 5.84829 + 2.53856 I, -8.59995 - 5.10712 I}, {0. + 0. I, 5.66679 + 16.9373 I, -6.42363 + 12.6433 I, -3.46984 - 4.99463 I, 9.86177 + 5.89218 I}, {0. + 0. I, 0. + 0. I, 12.6587 + 1.16017 10^-15 I, 1.67762 + 1.0062 I, 11.1778 + 4.4386 I}, {0. + 0. I, 0. + 0. I, 0. + 0. I, -13.136 + 5.90551 10^-16 I, -0.186487 + 0.0305276 I}, {0. + 0. I, 0. + 0. I, 0. + 0. I, 0. + 0. I, -7.85626 + 1.08747 10^-17 I} }
Results test out as correct (did not have a unitary matrix for the first example), but obvious variability. The failing doctest has a large condition number (~1000), but so does another test which is not failing. Notice that super-diagonal starts to agree at the right side.
It would be easy to fix doctests, by focusing on what the matrices should do, not on what they are. The eigenvalues could be stripped from the diagonal of the upper-triangular matrix to test as a list, since they should not vary.
This function seems too useful elsewhere (testing Hermitian, normal, etc) to abandon it. Thoughts on moving forward? We could
- Abandon failing test.
- Fix tests to concentrate on properties and not on uniqueness.
- Solicit testing on sage-devel to see how variable this is across systems and processors.
comment:4 Changed 10 years ago by
- Cc jason added
comment:5 Changed 10 years ago by
I wouldn't support abandoning the test. But of cause it tests something, that is not guaranteed by the (mathematical) definitions. My suggestion is: 1) Emphasize in the documentation that the decomposition cannot be assumed to be unique. In particular, it can even vary on different machines. 2) Actually the crucial: U is unitary, A is upper tridiagonal, the eigenvalues are in "right" order. 3) One can also test the Frobenius norm. That's not very deep, but it is a further invariant for the nilpotent part of A.
Your thoughts?
comment:6 in reply to: ↑ 3 Changed 10 years ago by
Replying to rbeezer:
- Fix tests to concentrate on properties and not on uniqueness.
+1
- Solicit testing on sage-devel to see how variable this is across systems and processors.
Also, you might write to the scipy list and ask how much variability we should expect on different platforms.
comment:7 Changed 10 years ago by
Hi Martin and Jason,
Thanks - those are all great suggestions. I'll probably get back to this tomorrow.
Rob
Changed 10 years ago by
comment:8 follow-up: ↓ 12 Changed 10 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
v2 patch adjusts the tests to not rely on exact values in the unitary matrix, nor above the diagonal in the upper-triangular matrix. It adds a few more checks on the diagonal entries, which are the eigenvalues of the matrix. There are also a couple of new checks on the 2-norm of the similar matrices, one for the real case and one for the complex case.
This passes tests on my desktop (Ubuntu, i7-2600) and on sage.math.
Jason - would you mind testing this on a Mac before we turn it loose?
Rob
comment:9 Changed 10 years ago by
- Status changed from needs_review to positive_review
That's very good. And now all tests pass.
I just noticed, and include this as a comment: Shouldn't we think of Rstyfying the matrix_double_dense file? Currently all your nice documentation doesn't show up in the reference manual. In particular, the schur method, that is somewhat unique to this class, won't be advertised.
comment:10 Changed 10 years ago by
Thanks, Martin! Yes, I went through this whole file once, but then started changing and adding whole routines and realized I was walking all over my formatting changes/updates. So I thought it would be best to wait just a little bit.
Once changes settle down I will go through and do some documentation (and minor code) clean-up and add this to the docs. It'll go up on #8046 once it happens. I've been checking my new documentation in the notebook as I go, so it should all format properly.
(I may add an exact Schur method, just for instructional use.)
comment:11 Changed 10 years ago by
- Merged in set to sage-4.7.alpha4
- Resolution set to fixed
- Status changed from positive_review to closed
comment:12 in reply to: ↑ 8 ; follow-up: ↓ 16 Changed 10 years ago by
Replying to rbeezer:
Jason - would you mind testing this on a Mac before we turn it loose?
It looks like this didn't happen, and I'm getting doctest failures. On a Mac (OS X 10.6.7, Intel):
File "/Applications/sage_builds/clean/sage-4.7.alpha4/devel/sage-main/sage/matrix/matrix_double_dense.pyx", line 1538: sage: T.round(4) Expected: [-13.5698 0.0 0.0 0.0] [ 0.0 -0.8508 -0.0 -0.0] [ 0.0 0.0 7.7664 0.0] [ 0.0 0.0 0.0 11.6542] Got: [-13.5698 0.0 0.0 0.0] [ 0.0 -0.8508 0.0 0.0] [ 0.0 0.0 7.7664 0.0] [ 0.0 0.0 0.0 11.6542]
On t2.math.washington.edu (Solaris on sparc):
File "/scratch/palmieri/clean/sage-4.7.alpha4/devel/sage/sage/matrix/matrix_double_dense.pyx", line 1538: sage: T.round(4) Expected: [-13.5698 0.0 0.0 0.0] [ 0.0 -0.8508 -0.0 -0.0] [ 0.0 0.0 7.7664 0.0] [ 0.0 0.0 0.0 11.6542] Got: [-13.5698 0.0 0.0 0.0] [ 0.0 -0.8508 0.0 -0.0] [ 0.0 0.0 7.7664 -0.0] [ 0.0 0.0 0.0 11.6542]
comment:13 Changed 10 years ago by
- Merged in sage-4.7.alpha4 deleted
- Resolution fixed deleted
- Status changed from closed to new
Unmerging this in sage-4.7.alpha5...
comment:14 Changed 10 years ago by
- Status changed from new to needs_work
comment:15 Changed 10 years ago by
On OpenSolaris 06/2009 (quad core Xeon):
File "/export/home/drkirkby/sage-4.7.alpha4/devel/sage/sage/matrix/matrix_double_dense.pyx", line 1538: sage: T.round(4) Expected: [-13.5698 0.0 0.0 0.0] [ 0.0 -0.8508 -0.0 -0.0] [ 0.0 0.0 7.7664 0.0] [ 0.0 0.0 0.0 11.6542] Got: [-13.5698 0.0 0.0 0.0] [ 0.0 -0.8508 0.0 -0.0] [ 0.0 0.0 7.7664 0.0] [ 0.0 0.0 0.0 11.6542] **********************************************************************
comment:16 in reply to: ↑ 12 ; follow-ups: ↓ 17 ↓ 19 Changed 10 years ago by
Replying to jhpalmieri:
Replying to rbeezer:
Jason - would you mind testing this on a Mac before we turn it loose?
It looks like this didn't happen, and I'm getting doctest failures. On a Mac (OS X 10.6.7, Intel):
Thanks, John and Dave.
I'm more concerned about #10837, so if you see those variable results, please let me know about that.
John - I have an account on Skynet, so I should do some testing there, I see. Do you have any pointers to etiquette/practices there? Do I build in my home directory? Should I "nice" things or limit the number of cores?
Rob
comment:17 in reply to: ↑ 16 ; follow-up: ↓ 18 Changed 10 years ago by
Some clarification of how this got into sage-4.7.alpha4: when I build a new Sage alpha release, it is tested on the buildbot which tests on many machines, including (most of) the Skynet machines. So this doctest problems did show up. However, in my original sage-4.7.alpha4, I also had #10837 and I mistook those doctest failures to come from #10837 instead. So I unmerged #10837 and went ahead to release sage-4.7.alpha4 without testing again. Experience shows that testing again after unmerging a patch from a candidate alpha is usually not necessary and only slows down things. In this case, it would have caught the errors from #11027.
Long story short: even if you don't test yourself, a patch will not get in a final Sage release if it fails doctests on Skynet.
comment:18 in reply to: ↑ 17 Changed 10 years ago by
Replying to jdemeyer:
Some clarification of how this got into sage-4.7.alpha4:
Thanks, Jeroen. No problems here.
I'd misunderstood and thought #10837 was in as well, so this helped clear that up for me.
I just need to get into the habit of testing on a wider variety of platforms if I'm going to tackle these numerical procedures. sage.math, of course, but I should add SkyNet
to my repertoire also.
Sorry for the trouble - I should have realized these -0.0 entries were a bad sign.
As always, thanks for all your careful work as release manager.
Rob
comment:19 in reply to: ↑ 16 Changed 10 years ago by
Replying to rbeezer:
John - I have an account on Skynet, so I should do some testing there, I see. Do you have any pointers to etiquette/practices there? Do I build in my home directory? Should I "nice" things or limit the number of cores?
Ignore the request, John. I think I have found my way. Maybe.
Rob
comment:20 Changed 10 years ago by
- Description modified (diff)
- Status changed from needs_work to needs_review
Well, it looks like this was just carelessness on my part, and not any kind of platform-variability problem. On other doctests, which seem to pass for the various testers, I have code like
sage: T = T.zero_at(1.0e-12).change_ring(RDF) sage: T.round(6)
The first line was missing on this one problem test - the doctest patch just adds in such a line. Likely a cut/paste when there should have been a copy/paste.
With the new patch, this passes tests on my 64-bit Ubuntu machine and on sage.math. I'm in the process of testing on SkyNet/cleo
. I'm confident this fix should be OK, but "once burned, twice shy."
This is ready for review, but I'd like to have the cleo result before it gets flipped to positive review.
Some discussion is necessary: As a doctest I get
Mathematica gives:
That is three different results, probably depending on the machine you run it on.