solve_right with matrix right hand side fails over RDF and CDF
fredrik-johansson opened this issue · 18 comments
Solving AX=B for a matrix B using solve_right or \ works over many rings, but not RDF or CDF. Example:
sage: MatrixSpace(QQ,10,10,).random_element() \ MatrixSpace(QQ,10,1).random_element()
[ -783/1145]
[ 15697/18320]
[ -4647/18320]
[ -2599/7328]
[ 20289/73280]
[ -3791/18320]
[-70107/36640]
[-12407/36640]
[ -1588/1145]
[ -7059/9160]
sage: MatrixSpace(RR,10,10,).random_element() \ MatrixSpace(RR,10,1).random_element()
[ 0.0620590489221758]
[ -0.584548090301576]
[ -0.106165379993850]
[ 1.52004291094317]
[ -1.03573096808637]
[ 0.0822404955372452]
[ -0.145002162865055]
[ 0.262055581969539]
[ 1.35298542346484]
[-0.0727293754429877]
sage: MatrixSpace(RDF,10,10,).random_element() \ MatrixSpace(RDF,10,1).random_element()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-16-0286e0222123> in <module>()
----> 1 MatrixSpace(RDF,Integer(10),Integer(10),).random_element() * BackslashOperator() * MatrixSpace(RDF,Integer(10),Integer(1)).random_element()
/home/fredrik/sage-6.3/local/lib/python2.7/site-packages/sage/misc/preparser.pyc in __mul__(self, right)
1413 (0.0, 0.5, 1.0, 1.5, 2.0)
1414 """
-> 1415 return self.left._backslash_(right)
1416
1417
/home/fredrik/sage-6.3/local/lib/python2.7/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix._backslash_ (build/cythonized/sage/matrix/matrix2.c:4193)()
/home/fredrik/sage-6.3/local/lib/python2.7/site-packages/sage/matrix/matrix_double_dense.so in sage.matrix.matrix_double_dense.Matrix_double_dense.solve_right (build/cythonized/sage/matrix/matrix_double_dense.c:11896)()
TypeError: vector of constants over Real Double Field incompatible with matrix over Real Double Field
Component: linear algebra
Author: Peter Wicks Stringfield, Jeroen Demeyer, Markus Wageringel
Branch/Commit: 7b07a27
Reviewer: Vincent Delecroix
Issue created by migration from https://trac.sagemath.org/ticket/17405
I attached a patch to this ticket.
Following the pattern set by the solve_right and solve_left in src/sage/matrix/matrix2.pyx, I tweaked the code in src/sage/matrix/matrix_double_dense.pyx to make solve_right and solve_left for RDF and CDF matrices accept matrices for b. If b is given as a matrix, the answer x will be returned as a matrix.
Please note:
- The solve_left in matrix2 uses solve_right, to reduce redundancy; while the solve_left in matrix_double_dense does NOT use solve_right, and instead duplicates some code.
- The solve_left and solve_right in matrix_double_dense are now more lenient than those in matrix2, since they accept not just vectors and matrices for b, but also things that can be coerced into vectors (like Python lists).
I didn't try to fix either of those issues.
Why the restriction that the right hand side only has one column?
Because it did not even occur to me that the right hand side could have more than one. Thanks for pointing that out.
I let b, the right hand side, have more columns.
I fixed a bug that is in apparently dead code and wouldn't matter even if it wasn't.
I made a branch for this ticket because that seems to be the proper way to do things here.
Also it seems like for non-square matrices the RDF/CDF right_solve and left_solve could just lean on the generic right_solve and left_solve? (The scipy solver needs square matrices, so right now the RDF/CDF code rejects them.) Anyway that would be for a different ticket if it even matters.
EDIT: I think I may have made the branch wrong? Sorry :(.
EDIT2: Figured it out. Branch is now correct.
Commit: 5dde1c0
Reviewer: Vincent Delecroix
Hello,
- You should not try to change the ring if it is already the good one as it performs a copy of the matrix
sage: m = matrix(ZZ,3,4,range(12))
sage: m.change_ring(ZZ) is m
False
sage: %timeit m.change_ring(ZZ)
1000000 loops, best of 3: 988 ns per loop
sage: m = matrix(ZZ,10,10,range(100))
sage: %timeit m.change_ring(ZZ)
100000 loops, best of 3: 1.51 µs per loop
This is bad since linear system might involve huge matrices.
- In your example, there is no need to print the input matrices
AandB(though it is not very bad).
Vincent
Having a look at this.
Author: Peter Wicks Stringfield, Jeroen Demeyer
Changed branch from u/peterwicksstringfield/17405_solve_right to u/jdemeyer/17405_solve_right
Replying to @videlec:
Hello,
- You should not try to change the ring if it is already the good one as it performs a copy of the matrix
I have updated the branch to use the coercion framework to figure out a common parent for A and B and made some minor improvements to the documentation.
I also added a deprecation for the case in which the argument b is not a vector or matrix, which was undocumented behavior and just complicates things for little gain.
In the future, solve_right for RDF/CDF should be updated to also handle the case of a non-square matrix A, which is already supported by generic matrices.
New commits:
bd6e4fc | Merge tag '9.1.beta4' into #17405 |
7b07a27 | 17405: fix solve_right and solve_left over RDF and CDF |
Changed branch from u/jdemeyer/17405_solve_right to u/gh-mwageringel/17405
Changed author from Peter Wicks Stringfield, Jeroen Demeyer to Peter Wicks Stringfield, Jeroen Demeyer, Markus Wageringel
Changed branch from u/gh-mwageringel/17405 to 7b07a27