Sage: Ticket #29226: Meta-ticket: issues with division of matrices, solve_left, solve_right
https://trac.sagemath.org/ticket/29226
<p>
This ticket collects several issues related to the division operation of matrices.
</p>
<ul><li><a class="closed ticket" href="https://trac.sagemath.org/ticket/29257" title="defect: use solve_left for division operation of matrices (closed: fixed)">#29257</a> <code>__truediv__</code> does not work for matrices of different dimensions:
<pre class="wiki">sage: set_random_seed(1)
sage: B, A = matrix.random(QQ, 2, 4), matrix.random(QQ, 4, 4)
sage: B * ~A; # works
sage: B / A # fails
...
TypeError: unsupported operand parent(s) for /: 'Full MatrixSpace of 2 by 4 dense matrices over Rational Field' an
d 'Full MatrixSpace of 4 by 4 dense matrices over Rational Field'
</pre>The corresponding <code>_backslash_</code> operation usually works.
</li></ul><ul><li><a class="closed ticket" href="https://trac.sagemath.org/ticket/29257" title="defect: use solve_left for division operation of matrices (closed: fixed)">#29257</a> For two square matrices, the default <code>__truediv__</code> implementation <code>left * ~right</code> computes the inverse of a matrix. This is unsuitable for numerical computations at least. A better implementation would invoke <code>solve_left</code> for matrices and vectors (<code>_backslash_</code> already calls <code>solve_right</code>).
</li></ul><ul><li><a class="closed ticket" href="https://trac.sagemath.org/ticket/12406" title="defect: solve_left and solve_right should use coercion (closed: fixed)">#12406</a> <code>solve_left</code> and <code>solve_right</code> should use coercion to find compatible parents
</li></ul><ul><li><a class="closed ticket" href="https://trac.sagemath.org/ticket/12406" title="defect: solve_left and solve_right should use coercion (closed: fixed)">#12406</a> <code>solve_right</code> and <code>solve_left</code> over <code>RDF</code>/<code>CDF</code> only work for square matrices, even though generic matrices support non-square matrices:
<pre class="wiki">sage: set_random_seed(0)
sage: A, B = matrix.random(QQ, 4, 3), matrix.random(QQ, 4, 2)
sage: A.solve_right(B); # works
sage: A.change_ring(RDF).solve_right(B.change_ring(RDF)) # fails
...
NotImplementedError: coefficient matrix of a system over RDF/CDF must be square, not 4 x 3
</pre>This could be solved by using <code>scipy.linalg.lstsq</code> to find a least-squares solution.
</li></ul><ul><li><a class="closed ticket" href="https://trac.sagemath.org/ticket/17405" title="defect: solve_right with matrix right hand side fails over RDF and CDF (closed: fixed)">#17405</a> <code>solve_right</code> and <code>solve_left</code> over <code>RDF</code>/<code>CDF</code> fail with matrix right hand side
</li></ul><ul><li><a class="closed ticket" href="https://trac.sagemath.org/ticket/13932" title="defect: solve_right fails with floating-point matrices (closed: duplicate)">#13932</a> <code>solve_right</code> and <code>solve_left</code> fail with floating point matrices (due to exact checking of the solution)
</li></ul><ul><li><a class="new ticket" href="https://trac.sagemath.org/ticket/28586" title="defect: broken solving a linear system over GF(3) (new)">#28586</a> <code>solve_right</code> and <code>solve_left</code> fail for <code>matrix_modn_sparse</code>
</li></ul><ul><li><a class="new ticket" href="https://trac.sagemath.org/ticket/29729" title="enhancement: improve solve_right for some inexact rings (new)">#29729</a> improve checking the result of <code>solve_right</code> for some inexact rings
</li></ul>en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/29226
Trac 1.1.6gh-mwageringelThu, 20 Feb 2020 19:46:01 GMT
https://trac.sagemath.org/ticket/29226#comment:1
https://trac.sagemath.org/ticket/29226#comment:1
<p>
What is a good solution for the first problem? Should there be inverse left/right actions of matrix multiplication, so that the coercion framework can be used to find these actions?
</p>
<p>
The implementation of <code>__truediv__</code> needs to be overwritten for matrices in any case in order to make use of <code>solve_left</code>. However, <code>solve_left</code> itself also needs to use some form of coercion to make sure the parents of the arguments are compatible (in <a class="closed ticket" href="https://trac.sagemath.org/ticket/17405" title="defect: solve_right with matrix right hand side fails over RDF and CDF (closed: fixed)">#17405</a>, I implemented this in an ad-hoc way which could be copied). I am not sure whether implementing inverse actions would help with this or make things needlessly complicated.
</p>
Ticketgh-mwageringelFri, 28 Feb 2020 13:17:59 GMTdescription changed
https://trac.sagemath.org/ticket/29226#comment:2
https://trac.sagemath.org/ticket/29226#comment:2
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/29226?action=diff&version=2">diff</a>)
</li>
</ul>
<p>
The coercion is now implemented in <a class="closed ticket" href="https://trac.sagemath.org/ticket/12406" title="defect: solve_left and solve_right should use coercion (closed: fixed)">#12406</a>, without the need for inverse actions. The <code>__truediv__</code> operation is implemented in <a class="closed ticket" href="https://trac.sagemath.org/ticket/29257" title="defect: use solve_left for division operation of matrices (closed: fixed)">#29257</a>.
</p>
TicketmjoSun, 05 Apr 2020 21:49:59 GMTcc set
https://trac.sagemath.org/ticket/29226#comment:3
https://trac.sagemath.org/ticket/29226#comment:3
<ul>
<li><strong>cc</strong>
<em>mjo</em> added
</li>
</ul>
<p>
While technically we can define "division" to mean right-multiplication-by-inverse in a noncommutative ring in a mathematically consistent way, I'm not sure we should have -- particularly for matrices. I have never written <code>A/B</code> with matrices except by mistake, and I've never seen anyone else do it on purpose either. I would much rather get an error if I ever do it again (<em>especially</em> if their dimensions are mismatched) than silently get a somehow-correct-but-definitely-not-what-I-meant answer back. In other words, I think we should just kill matrix-matrix truediv. Having a pedantically correct definition that does the wrong thing 100% of the time isn't helping anything.
</p>
Ticketgh-mwageringelMon, 06 Apr 2020 19:40:53 GMT
https://trac.sagemath.org/ticket/29226#comment:4
https://trac.sagemath.org/ticket/29226#comment:4
<p>
Well, the division operator for matrices has existed in Sage for years, so removing it at this point would be quite unfortunate. In <a class="closed ticket" href="https://trac.sagemath.org/ticket/29257" title="defect: use solve_left for division operation of matrices (closed: fixed)">#29257</a>, I am just trying to improve the implementation of this matrix operator, to make it as good as it can be.
</p>
<p>
Let me explain why this operator is important for matrices. Here is an example of an actual computation I did a few days ago
</p>
<pre class="wiki">S \ U.H * A * U / S
</pre><p>
where <code>U,S,A</code> are matrices. In Sage, this currently only works if <code>U</code> is square, which is unnecessarily restrictive, so you are forced to write
</p>
<pre class="wiki">~S * U.H * A * U * ~S
</pre><p>
This comes with its own problems because it is slower and numerically less accurate due to the inversion of a matrix. A better way to compute this would be
</p>
<pre class="wiki">S.solve_left(S.solve_right(U.H) * A * U)
</pre><p>
Writing it like this is rather backwards and difficult to remember. After <a class="closed ticket" href="https://trac.sagemath.org/ticket/29257" title="defect: use solve_left for division operation of matrices (closed: fixed)">#29257</a>, this will be equivalent to the first example, which is a much easier syntax as the binary operators compose nicely.
</p>
<p>
The slash and backslash operators are very standard notation in Matlab, so many people will be familiar with it. After all, it is Sage's mission to be a viable alternative to Matlab among others. Considering all our work on <a class="closed ticket" href="https://trac.sagemath.org/ticket/12406" title="defect: solve_left and solve_right should use coercion (closed: fixed)">#12406</a> to improve the backslash operator, I am a bit surprised by your concerns about the forward slash operator. I have difficulties to see how this is doing the wrong thing 100% of the time.
</p>
<p>
As for the numerical issues, here is an example from the <a class="ext-link" href="https://www.mathworks.com/help/matlab/ref/inv.html#bu6sfy8-1"><span class="icon"></span>Matlab documentation</a> that shows that using the slash operator to solve <code>A*x = b</code> is typically faster and more accurate than computing the matrix inverse:
</p>
<pre class="wiki">sage: import numpy
sage: n = 500
sage: Q = matrix.random(RDF, n).QR()[0]
sage: d = numpy.logspace(0, -10, n)
sage: x = matrix.random(RDF, n, 1)
sage: A = Q * matrix.diagonal(CDF, d) * Q.H
sage: b = A * x
sage: y = ~A * b # first solution
sage: z = A \ b # second solution
sage: %timeit y = ~A * b # slower
10 loops, best of 5: 29.6 ms per loop
sage: %timeit z = A \ b # faster
100 loops, best of 5: 13.7 ms per loop
sage: (A * y - b).norm() # large error
1.2981997570563244e-07
sage: (A * z - b).norm() # small error
1.7440675515723775e-15
</pre>
TicketmjoMon, 06 Apr 2020 23:13:57 GMT
https://trac.sagemath.org/ticket/29226#comment:5
https://trac.sagemath.org/ticket/29226#comment:5
<p>
I suspect it's a religious matter, don't take me too seriously. Of your two examples,
</p>
<pre class="wiki">S \ U.H * A * U / S
</pre><pre class="wiki">S.solve_left(S.solve_right(U.H) * A * U)
</pre><p>
I would greatly prefer the second, because I can tell what it does. The first depends on me knowing what matrix division does, how it associates, and how it's implemented. Whereas the second example is clearly just solving two systems, and I can tell you which two systems those are right away.
</p>
<p>
What it comes down to for me is that, unless you're a heavy MATLAB user, division of matrices is basically a nonsense operation. I would never think to divide two matrices, so if <code>A / B</code> appears in my own code, it's a mistake. And I'm just guessing based on personal experience (again, to be taken with a grain of salt) that most users of Sage are the same.
</p>
<p>
This is reminiscent of the <a class="missing wiki">Foldable/Travserable?</a> proposal in Haskell. If you open up a Haskell prompt and type <code>length (1,2)</code>, you'll get back <code>1</code> as the answer. Why? Becuase if you invest 20 years of education and ultimately wind up studying category theory, you'll see that there's a perfectly canonical way to define "fold" on the category of tuples
and if you implement <code>length</code> in a perfectly canonical way, well, <code>1</code> is the only possible length that <code>(1,2)</code> could have! This is of course ridiculous: the job of the compiler is to prevent me from writing <code>length (1,2)</code> when I mean <code>length [1,2]</code>, and enabling <code>length</code> to work on a data type that it's not going to act sanely on prevents it from doing so. I see matrix division the same way.
</p>
<p>
The fact that MATLAB implements it is a good argument, but MATLAB code is also impossible to debug for this same reason. If my algorithm has a bug and tries to multiply a PNG image by a matrix of strings, MATLAB might silently return a column vector of TCP/IP packets. Tiny errors thus propagate, and make it impossible to localize the cause of a bug.
</p>
<p>
Sage (via python) isn't even statically typed, so it's easier to make the argument here than in Haskell that if the data type supports <code>__truediv__</code>, then we should use it. It's also very possible that everyone else just plain disagrees with me, and that's fine too. I'm just throwing in my two cents: it's easier for me to write big programs if the compiler keeps me from doing things I don't mean to do. And for me, dividing by a matrix is something I never mean to do.
</p>
TicketmkoeppeThu, 07 May 2020 02:29:39 GMTmilestone changed
https://trac.sagemath.org/ticket/29226#comment:6
https://trac.sagemath.org/ticket/29226#comment:6
<ul>
<li><strong>milestone</strong>
changed from <em>sage-9.1</em> to <em>sage-9.2</em>
</li>
</ul>
Ticketgh-mwageringelSat, 23 May 2020 19:01:17 GMTdescription changed
https://trac.sagemath.org/ticket/29226#comment:7
https://trac.sagemath.org/ticket/29226#comment:7
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/29226?action=diff&version=7">diff</a>)
</li>
</ul>
Ticketgh-kliemSat, 15 Aug 2020 10:46:38 GMT
https://trac.sagemath.org/ticket/29226#comment:8
https://trac.sagemath.org/ticket/29226#comment:8
<p>
@mjo: I positively reviewed <a class="closed ticket" href="https://trac.sagemath.org/ticket/29257" title="defect: use solve_left for division operation of matrices (closed: fixed)">#29257</a>, because I think it is an improvement in making right and left division behave consistently. The use of left division as an alias for solve right is even advertised in <a class="ext-link" href="https://doc.sagemath.org/html/en/tutorial/tour_linalg.html"><span class="icon"></span>https://doc.sagemath.org/html/en/tutorial/tour_linalg.html</a>.
</p>
<p>
Now, making them behave consistently can also mean that they are both undefined or both raise a warning. Maybe a longterm warning as "it is assumed that you meant to use solve left; it is recommended to use <code>solve_left</code>" would work, because it helps users to find the correct things coming from MATLAB. Whoever wants to, can just use it, but we would consider it unclean and not accept it e.g. for sage's code.
</p>
TicketmjoMon, 17 Aug 2020 12:59:24 GMT
https://trac.sagemath.org/ticket/29226#comment:9
https://trac.sagemath.org/ticket/29226#comment:9
<p>
(I don't really expect to win this argument.)
</p>
<p>
Having left- and right-division act consistently does makes sense if you think of the backslash operator as division-with-its-arguments-reversed. But backslash has two killer features: it already doesn't work in library code, and you can't use it to divide scalars, so there's no risk of concealing a type error.
</p>
TicketmkoeppeSat, 24 Oct 2020 20:15:01 GMTmilestone changed
https://trac.sagemath.org/ticket/29226#comment:10
https://trac.sagemath.org/ticket/29226#comment:10
<ul>
<li><strong>milestone</strong>
changed from <em>sage-9.2</em> to <em>sage-9.3</em>
</li>
</ul>
TicketmkoeppeThu, 15 Apr 2021 02:16:37 GMTmilestone changed
https://trac.sagemath.org/ticket/29226#comment:11
https://trac.sagemath.org/ticket/29226#comment:11
<ul>
<li><strong>milestone</strong>
changed from <em>sage-9.3</em> to <em>sage-9.4</em>
</li>
</ul>
TicketmkoeppeTue, 10 Aug 2021 17:30:20 GMTmilestone changed
https://trac.sagemath.org/ticket/29226#comment:12
https://trac.sagemath.org/ticket/29226#comment:12
<ul>
<li><strong>milestone</strong>
changed from <em>sage-9.4</em> to <em>sage-9.5</em>
</li>
</ul>
TicketmkoeppeSat, 18 Dec 2021 19:41:23 GMTmilestone changed
https://trac.sagemath.org/ticket/29226#comment:13
https://trac.sagemath.org/ticket/29226#comment:13
<ul>
<li><strong>milestone</strong>
changed from <em>sage-9.5</em> to <em>sage-9.6</em>
</li>
</ul>
Ticket