@double/adjoint: Add new function. #1192
Conversation
inst/@double/adjoint.m
Outdated
|
|
||
|
|
||
| function A = adjoint (M) | ||
| A = double (adjoint (vpa (M))); |
There was a problem hiding this comment.
If user sets digits(3) then this code will respect that and (I guess?) give a poor result even though they would not expect a double function to care...
Similarly, if user sets digits(1000) this will take a lot longer before truncating to double precision.
There was a problem hiding this comment.
Good catch! It's probably better to call Let's use pycall_sympy__ directly instead of relying on vpa.vpa (M, 16) instead.
inst/@double/adjoint.m
Outdated
| %! M = randn (2); | ||
| %! A = [M(2,2) -M(1,2); -M(2,1) M(1,1)]; | ||
| %! assert (isequal (adjoint (M), A)); | ||
| %! assert (isequal (adjoint (adjoint (M)), M)); |
There was a problem hiding this comment.
these look like floating point equality tests. Is it reasonable to expect a double-precision function to satisfy this test?
There was a problem hiding this comment.
Indeed, should I do abs (x - y) < 1e-15 instead?
There was a problem hiding this comment.
I suppose so, in a relative sense would be better. Maybe some of the other @double/<fcn> have some clues.
I worry a bit about the numerical stability of this sort of implementation (SymPy is not well-known for stability of its implementations when sympy.Float is involved: basically unless it calls mpmath I don't really trust it).
The short-term fix is that I think a random test is ill-advised unless you know the algorithm used is numerically stable. I'd suggest a couple hardcoded cases.
There was a problem hiding this comment.
How about a nearly singular one? [0 1; 0.0001 1]. Does that still give 1e-15?
There was a problem hiding this comment.
It works for this particular case in sympy 1.10.1 with the default method of computing determinant. But let's keep it simple and use a fixed constant to avoid any future issues...
Fixes gnu-octave#1049. * inst/@double/adjoint.m: New file. * INDEX: Adjust accordingly. * NEWS: Adjust accordingly.
Fixes #1049.