How Does mldivide Work for Full, Square, Full Rank, Triangular Matrix?

How does mldivide work when the input matrix A is full, square, full rank, and triangular? The doc page just say it uses a triangular solver, and I'd like more insight into what that triangular solver is.

An iterative approach can be used, e.g., for a lower triangular A

A = [1 0 0;2 3 0;4 5 6];
b = [1; 2; 3];
% A*x = b
% solve for x(i) iteratively, could be easily done with a for loop
x(1) = b(1)/A(1,1);
x(2) = (b(2) - A(2,1)*x(1)) / A(2,2);
x(3) = (b(3) - A(3,1)*x(1) - A(3,2)*x(2)) / A(3,3);
A*x(:) - b

ans = 3×1

0 0 0

Does mldivide use a solver that's more efficient and/or more numerically robust?

[Bruno Luong]

I'm not aware of any drawback of the standard forward/backward substitution algorithm, and any better method.

[John D'Errico]

Why should there be something more numerically robust? Whatever that means here. Just because you want it, this is not suffcient for something to exist. If you are forced to divide by zero, or some very small number, nothing in mathemagics will prevent that divide.

As for the iterative approach you write, that is great for a 3x3 matrix. For larger matrices, you will still be wanting to use linear algebraic operations to perform the computations, thus using the BLAS for efficiency.

Essentially, the code is doing what you wrote, but by use of the BLAS.

[Paul]

@Bruno Luong:

Thanks for the straightforward response.

@John D'Errico

Your comment makes inferences about my positions on what "should there be" and what I "want," neither of which are true. Perhaps you should edit your comment.

As to the substance of your comment, good idea to check the BLAS. Perhaps Matlab uses DTRSV for the double, real case (and STRSV and CTRSV for single and complex), though I don't see "linear algebraic operations" in DTRSV other than multiplication, subtraction, and division of scalars, which are the same operations in my admittedly simple-minded iterative approach.

[Steven Lord]

The doc page just say it uses a triangular solver, and I'd like more insight into what that triangular solver is.

Might I ask how you're hoping to use that information if it were available? How would your code be different if we told you we used approach X instead of approach Y?

I tend to discourage users from tightly coupling their code to the implementation the release they're using has, as it may mean that when or if we introduce approach Z that is generally "better" their code (tied to approaches X or Y) may not be able to take advantage of that improvement. Sure, some coupling is unavoidable or perhaps even good, but in this case I'd say the coupling is "Use backslash instead of trying to build your own solver."

[Paul]

Hi Steven,

I wouldn't use that information in the least. I'd use mldivide, \ without giving it a second thought for a triangular A should I ever have the need.

I agree with your comment to "use backslash instead of trying to build your own solver" for my work, and it's probably very good advice for most, if not all, people's production code.

However, I've seen such advice (sometimes more as an admonishment, though not from you in my recollection) given on this forum to people who might only be trying to learn about something, e.g., "Why are you bothering to write your own solver? Just the use the code that the experts provide!" One of the best ways to learn is to try to solve a problem oneself, IMO. And in these cases the native Matlab solution can serve as a tool for verification and comparison, which makes it a wonderful resource.

In this instance, I came across a problem in linear sytems that reduced to a triangular A*x = b problem where the author said such problems always have a solution, as one can see from the iterative substitution approach. But the point of the author was just to show that a solution always exists to the linear equations as opposed to showing how to best solve them. So I consulted mldivide, \ to see how Matlab would solve the triangular problem, which led to my question. So more out of curiosity than anything else.

[Bruno Luong]

@Steven Lord Sometime user needs a precise reference of the algorithm to analyze the the complexity, the performance, the memory required, etc... within a context of watever he and she is developping; especially if it is under scientific review.

[Bruno Luong]

@Paul "So I consulted mldivide, \ to see how Matlab would solve the triangular problem"

This is more related to how MATLAB compute the rank and disable the projection of the solution into what it considers to the kernel of the system.

In is sense the question is interesting for triangular system since I actually don't know how it estimates the rank.

For the non-triangular full matrix it performs QR decomposition with columns permutation, so the diagonal is decreses in absolute and it cut with the criteria I have asked in another recent post about rank revealing.

If the matrix is already triangular, actually I don't know how the rank and the kernel is actually estimated, since no QR would be performed on top of the triangular matrix, I believe the rank estimation must be performed before the dtrsv (or whatever the equivalent blas function MATLAB uses.

Look for a piece of algorithm won't allow to answer your question IMO.

[Paul]

Hi Bruno,

The comments in DTRSV states: "No test for singularity or near-singularity is included in this routine. Such tests must be performed before calling this routine." As you suspected.

I was wondering about the rank aspect as well, but since it wasn't central to what I was curious about I just specified full rank in my question to eliminate that aspect of the problem.

The flow chart on mldivide, \ just has the test "Is A triangular?" Are you thinking that that test would be more accurately stated as "Is A triangular to within some rank tolerance?"

For a triangular matrix, is the mathematical criterion for full rank all elements of the diagonal are nonzero?

[Bruno Luong]

I think the test of the matrix (triangular) form does not have tolerance > 0 (such as MATLAB test for symmetric or hermitian). However the rank testing must have ome sort of tolerance.

I have some doubt that MATLAB uses for strict mathematical criteria of full rank for practical numerical rank (for triangular matrix), that would be dangeruous. But this is just my assumption.

[Bruno Luong]

This code shows that MATLAB (R2023B) rank checking are not identical for triangular and not triangular matrix

small = eps(1);
T=[small 0; 
   1 1];
b = [1;0];
Tf = fliplr(T); % not a upper or lower triagular matrix

% QR method
x1=Tf\b % trigger warning

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.110223e-16.

x1 = 2×1

1.0e+15 * -4.5036 4.5036

x1=flipud(x1) % restore the order due to matrix column swap

x1 = 2×1

1.0e+15 * 4.5036 -4.5036

% back substitution (?)
x2=T\b % no rank warning for lower triangular matrix

x2 = 2×1

1.0e+15 * 4.5036 -4.5036

% results however matche
x1-x2

ans = 2×1

0 0

[Paul]

small = eps(1);
T=[small 0; 
   1 1];
b = [1;0];

No warning, even though rcond(T) < eps (I thought 0 < rcond(T) < eps is the criterion to trigger the warning)

format long e
rcond(T)

ans = 
     1.110223024625156e-16

rcond(T) < eps

ans = logical
   1

x1 = T\b

x1 = 2x1

1.0e+00 * 4.503599627370496e+15 -4.503599627370496e+15

Make T(1,1) a very small number

T(1,1) = T(1,1)/1e10

T = 2x2

2.220446049250313e-26 0 1.000000000000000e+00 1.000000000000000e+00

rcond(T)

ans = 
     1.110223024625157e-26

x2 = T\b

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.220446e-26.

x2 = 2x1

1.0e+00 * 4.503599627370496e+25 -4.503599627370496e+25

Does the existence of the warning imply that x2 is computed with a different solver than x1?

Also, why is RCOND in the warning message 2x larger than rcond(T)?

[Bruno Luong]

My guess is that the warning during solving x=A\b is not alway triggered based on matrix conditioning number or matrix norm. Those are not cheap to compute.

Rather it's based first on eventual intermediate quantity computed while the linear solver algorithm does the work: smallest and largest diagonal elements of A if A is triangular form, and diagonal of R if QR decomposition is used (generic full matrix).

This code seems to indicate that

1. substitution is used directly on T (small is compared to 1 for singularity detection) and
2. QR is used on Tf (small is compared to norm([1; 1]) = sqrt(2) for singularity detection)

[Christine Tobler]

Yes, rcond is not computed in x = A\b if A is triangular. This is because computing the approximate rcond value would be several times more expensive than solving the linear system itself. Instead, a heuristic based on the diagonal values of A is used.