ans =
Should mldivide Return a Solution for Square, Rank-Deficient, but Consistent Set of Linear Equations?
Ältere Kommentare anzeigen
load lindata
A and b define a linear system of equations Ax = b.
A is 12 x 12 with rank
rank(A)
But the linear equations are consistent with an infinite number of solutions. One such solution is
x1 = repmat([0;0;1],4,1);
Another solution can be found with the sym version of mldivide
x2 = A\b;
Verify
A*[x1,x2] - b
But the numeric version of mldivide fails
x3 = double(A)\double(b)
"
- If A is a square n-by-n matrix and B is a matrix with n rows, then x = A\B is a solution to the equation A*x = B, if it exists."
My interpretation of that last clause is that it would be properly read as "if a solution exists."
Here, many solutions exist but mldivide didn't find a solution. Is this an oversight on the doc page, or am I misreading it, or would it be reasonable to expect mldivide to return a solution for this sort of problem?
Akzeptierte Antwort
Christine Tobler
am 10 Mär. 2025
Hi Paul!
You can use
x = linsolve(A, b, struct('RECT', true));
to use the QR method of mldivide even if A is square.
To get the minimum-norm solution of x, instead of just any solution x like mldivide does, you can use
x = lsqminnorm(A, b);
To see what algorithm is chosen by mldivide, you can call
dA = decomposition(A);
and check the Type property of the returned object.
3 Kommentare
The QR approach in linsolve and mldivide (they are the same) is to do a permuted QR decomposition, estimate the rank r of A based on the diagonal of the R matrix, and then truncate the system to only use the first r rows and columns of R. See the Answer here for details.
For linsolve, setting RECT will have A treated like a rectangular matrix (meaning QR decomposition is used), which safe for a square matrix. If a non-upper triangular matrix is passed in, but A is not upper triangular, its lower triangle will just be assumed to be zero. The same is true for lower triangular and symmetric, where we just read one triangle of the matrix and assume the other is its conjugate transpose.
So the output is always "valid" in the sense that the operation is valid for how the input matrix is interpreted, but it could be confusing. I will let the documentation team know your suggestions.
For the statement in mldivide, perhaps that second bullet point should say "If A is a square n-by-n matrix and B is a matrix with n rows, then x = A\B is a solution to the equation A*x = B, if A is well-conditioned.".
On having mldivide detect this case: Performance is part of the problem, particularly because a majority of calls to mldivide are from applications that always use square matrices. But more importantly, the QR algorithm we use relies on a tolerance for detecting the rank of a matrix, which is fine for a least-squares problem, but may not give the expected result otherwise.
For example, if A was a diagonal matrix with very different scaling on the diagonal, the smallest elements on the diagonal would be cut off, with those elements of the output getting set to zero:
format shortg
linsolve(diag(10.^(0:-5:-30)), ones(7, 1), struct('RECT', true))
linsolve(diag(10.^(0:-5:-30)), ones(7, 1))
This is because the QR-based algorithm assumes that any change of A that is less than eps * norm(A) (more-or-less) is irrelevant, so these small diagonal elements should have no impact and should be removed to avoid potential Inf/NaN/overflow. But if we're not looking for a least-squares solution, but just for a direct solution, these elements have massive impact, because the smaller the element is, the larger its inverse.
Paul
am 11 Mär. 2025
Weitere Antworten (1)
Jack
am 8 Mär. 2025
0 Stimmen
Hey Paul,
Good question. mldivide (\) is supposed to find a solution when one exists, but when A is square and singular, MATLAB treats it as numerically unstable and returns NaNs instead of picking one solution.
The reason A\b works in sym but not in double is that the symbolic version does exact algebra, so it can return a valid solution even when there are infinitely many. The numeric version, on the other hand, works with floating-point numbers. Since A is singular, MATLAB can’t get a stable least-squares solution, so it just throws out NaNs.
If you still want a solution in double, you can use the pseudoinverse instead:
x3 = pinv(A) * b;
That should give you one possible solution.
As for the doc page, technically mldivide is doing what it says—it guarantees a solution when A is invertible or full-rank, but when it’s singular, it just gives a warning and returns NaNs instead of trying to pick one from infinitely many.
Follow me so you can message me anytime with future MATLAB questions. If this helps, please accept the answer as well.
13 Kommentare
Paul
am 8 Mär. 2025
Jack
am 8 Mär. 2025
mldivide (the backslash operator) is designed to return a solution only when it can compute a stable, unique solution. When A is square but rank‐deficient, the system has infinitely many solutions, and the numeric algorithm used in double precision cannot reliably pick one without risking numerical instability. In such cases, MATLAB deliberately returns NaNs rather than attempting to choose one solution arbitrarily. By contrast, the symbolic version performs exact algebra and can return one of the many possible solutions. If you need a solution in double precision for a consistent, but rank‐deficient system, you can use the pseudoinverse (pinv) as a workaround (for example, x = pinv(A)*b).
David Goodmanson
am 9 Mär. 2025
Bearbeitet: David Goodmanson
am 12 Mär. 2025
Hi Paul,
consider the following situation for Ax = b
% A (mxn) % x (nxq) % y (mxq)
m = 5; % m = 2,3, ...
n = 6;
q = 3;
A = rand(m,n);
A(end,:) = A(end-1,:) % make A rank deficient by 1
% A(:,end) = A(:,end-1); % another way to do it
b = rand(m,q);
x = A\b
Warning: Rank deficient, rank = 4, tol = 2.371612e-15.
x =
0.2542 0.6945 0.0126
-1.0880 -0.7025 -0.3020
0.0745 0.2533 0.4463
0 0 0
0.9049 0.4365 0.5911
0 0 0
This works for any number of rows, so that A can be wider than tall or taller than wide, with one exception. When m = 6 and A is square, you get the nan deal:
Warning: Matrix is singular to working precision.
x1 =
NaN NaN NaN
NaN NaN NaN
NaN NaN NaN
NaN NaN NaN
NaN NaN NaN
NaN NaN NaN
I guess the considerations are different when A is square, but I think it's interesting that backslash declines to come up with an answer in the square case but is fine doing so when A is rectangular.
Jack
am 9 Mär. 2025
mldivide follows different paths based on the shape of A, and when A is square, it defaults to LU, which breaks down for singular matrices. When A is rectangular, QR comes into play, which is why it sometimes works in rank-deficient cases.
Trying qr manually won’t help much since QR alone doesn’t solve the system—it just decomposes A. But your idea of mldivide checking for rank and falling back to pinv(A) * b when the system is consistent could work in some cases. The challenge is handling A*x = B when some columns of B are consistent and others aren’t. MATLAB would have to decide on a per-column basis, which might make the implementation more complex.
I agree that the doc could be clearer—saying "usually a solution" instead of an absolute statement would better reflect edge cases like this. And having an option like [x, alg] = mldivide(A, b) to return the solver used would be useful, especially when dealing with singular or ill-conditioned matrices.
John D'Errico
am 9 Mär. 2025
Bearbeitet: John D'Errico
am 9 Mär. 2025
I agree with most of what you said. But I disagree that MLDIVIDE should check for rank, and if that suggests a singular system, should default to the PINV solution. The first reason I see is PINV will sometimes be highly computationally intensive. And as you point out, if some columns of b live in the column space of A (and so are consistent) then the use of PINV could be confusing. Even worse is the case where a column of b lives almost in the column space of A.
A = magic(4) % rank = 3
b = A(:,1)
A\b % backslash actually gets it right, but pinv will certainly give some other result
pinv(A)*b
b2 = [A(:,1), A(:,1) + 1e-15*null(A)]
The columns of B are almost identical, but the second column lies just slightly out of the column space of A. Should MLDIVIDE decide to use PINV on the second column, but not the first? UGH.
I very much agree the documentation could be made more clear.
Walter Roberson
am 9 Mär. 2025
Bearbeitet: Walter Roberson
am 9 Mär. 2025
I suspect you wanted to proceed on to a demonstration with b2 ?
A = magic(4) % rank = 3
b2 = [A(:,1), A(:,1) + 1e-15*null(A)]
A\b2
pinv(A)*b2
John D'Errico
am 9 Mär. 2025
Oh, yes. Forgot that, :)
Anyway, I honestly think the real answer should be that when backslash sees a singular matrix, if should issue a fluorescent warning, of some different color. It should say something to the effect of DON'T TRUST ANYTHING YOU SEE FROM THIS COMPUTATION!
Jack
am 10 Mär. 2025
Pinv becomes computationally intensive when dealing with large or nearly singular matrices because it relies on the singular value decomposition (SVD), which has a significantly higher computational cost compared to the LU or QR decompositions used in mldivide for full-rank matrices. In cases where A is large or ill-conditioned, the SVD (and thus pinv) can be very slow.
Regarding the idea of adding a second output to mldivide that indicates which algorithm was used (e.g., LU, QR, or a fallback to pinv), it's an interesting suggestion for debugging and transparency. However, mldivide is designed to be a simple, one-line operator that hides its internal decision-making process. Changing its interface to return diagnostic information would complicate its use and could introduce backward compatibility issues. A better approach might be to offer a separate diagnostic function that reports on the condition of A and details the decomposition strategy chosen, without altering the behavior of mldivide itself.
John D'Errico
am 10 Mär. 2025
Bearbeitet: John D'Errico
am 10 Mär. 2025
Under what circumstances is pinv computationally intensive? Large arrays.
A = rand(2000);timeit(@() lu(A))
ans =
0.026602
A = rand(2000);timeit(@() svd(A,0))
ans =
0.84203
2Kx2K is not even what I would call large in this context these days.
Would users freak out, if mldivide started to take 30 or 40 times as much time to return a result? On a problem where the result is not even unique?
Instead, just return a warning, that says, IF they understand what the warning means, that this result is suspect.
Should MLDIVIDE return the specific algorithm used, IF asked? I think more information is never a bad thing. And only someone who understands what that information means will probably ever ask that "question", at least in this case. The only issue then is how they would indicate that information. A simple algorithmic name {"QR", "LU", CHOL"} might be insufficient to some, but they could always read the help docs to expand the meaning.
Pinv relies on computing the singular value decomposition (SVD), which has cubic complexity and can be significantly slower than LU or QR for large matrices. In a 2000×2000 matrix, LU takes roughly 0.03 seconds while SVD takes about 0.84 seconds—about 30 times slower. If mldivide were modified to return a solution via pinv for square, singular matrices, users might indeed notice a substantial performance penalty, especially on larger problems.
In essence, mldivide is designed to return a solution only when it can do so reliably using its standard LU (or QR, for non-square cases) methods. When the matrix is singular and the system has infinitely many solutions, returning a particular solution via pinv (which computes the minimum-norm solution) could be computationally intensive and might not be what a typical user expects, given that the solution isn’t unique.
Regarding returning diagnostic information about which algorithm was used, that could be useful for advanced users. However, mldivide’s simplicity as a one-line operator is one of its strengths, and changing its interface (for example, adding a second output to indicate the algorithm) could break backward compatibility or complicate usage for the majority of users. A more balanced approach might be to provide a separate diagnostic function that inspects the matrix condition and details the decomposition strategy used internally without altering mldivide’s behavior.
Walter Roberson
am 12 Mär. 2025
The vast majority of mldivide use is through the \ operator, which is a syntax that does not permit returning multiple outputs. In order for it to be used in a syntax that permitted returning multiple outputs, it would pretty much have to be invoked as
[Result, ExtraInfo] = mldivide(A,b);
though I guess it would also work to invoke
[Result, ExtraInfo] = A\b;
Currently that syntax generates an error about "Too many output arguments."
So there cannot be any present use of returning multiple outputs from \ or mldivide() that does not generate an error message. Accordingly, the only possible interferance with backwards compatibility by adding a second output argument would be against people deliberately writing
try; [Result, Extrainfo] = A\b; report_that_forbidden_event_has_occurred(); end
the number of which is surely vanishingly small.
Kategorien
Mehr zu Linear Algebra finden Sie in Hilfe-Center und File Exchange
am 12 Mär. 2025
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
