I've been tinkering a little in Fortran (2008) and wrote the following to solve $Rx=b$ for $R\in\mathbb{R}^{n\times n}$ upper-triangular, $x,b\in\mathbb{R}^n$.

My code looks like this:

function backsub_upper(R, b) result(x)

    ! R is a square upper triangular matrix
    ! b is the vector to be solved for
    ! x is the solution to Rx=b

    real(kind=wp), dimension(:,:), intent(in)   :: R
    real(kind=wp), dimension(:) , intent(in)    :: b

    real(kind=wp), dimension(:)                 :: x(size(b)), bc(size(b))
    integer(kind=sp)                            :: i,j

    x = 0.0_wp
    bc(:) = b(:)
    do j=1,size(bc)
        i = size(bc)+1-j
        x(i) = bc(i)/R(i,i)

        if (i/=1) bc(1:i-1) = bc(1:i-1) - x(i)*R(1:i-1,i)
    end do

end function backsub_upper

I tested with random $R$ (upper triangular) and random $b$. Works great but for more than n=~50 the error starts going crazy. I could keep the error down by enforcing some diagonal dominance, but just wondering if there is a better algorithm that control the error better for general matrices?

I realise random matrices are not 'general' but this feature got me worried

  • $\begingroup$ What did you compare your computed solution with to assess its accuracy, if you took random $R$ and $b$? $\endgroup$ May 21 '20 at 17:15
  • $\begingroup$ I checked the norm2(matmul(R,x)-b) $\endgroup$ May 23 '20 at 10:08
  • $\begingroup$ This is a bit difficult because actually I am implementing a back substitution of course to do solve upper triangular matrices such as generated from QR decomposition. So @rchilton1980 answer is meaningful. On the other hand Federico Poloni's answer tells me "No. This is as good as it gets". I was actually wondering if I had implemented this in a way that exaggerates round-off errors. Not sure which answer to accept here. $\endgroup$ May 23 '20 at 10:13
  • 1
    $\begingroup$ That is usually called residual, not error. What do you mean exactly by 'going crazy'? Some growth in the residual with the dimension is to be expected, because $R$ and $b$ have more and larger entries. The recommended thing to check is the relative residual norm2(matmul(R,x)-b) / norm2(R)norm2(x), or also norm2(matmul(R,x)-b) / norm2(b) which should be of the same order unless $R$ and $b$ are chosen in a special way. This quantity should always be in the ballpark of machine precision (times a factor $n$ or so), irrespective of condition numbers. $\endgroup$ May 23 '20 at 11:53
  • $\begingroup$ netlib.org/lapack/lapack-3.1.1/html/dtrsv.f.html $\endgroup$
    – vibe
    May 23 '20 at 12:23

The issue here is probably the choice of input, a random triangular matrix can have really bad conditioning. A better choice here is the L or U factor from the LU decomposition of a random matrix, or the R factor from the QR decomposition of a random matrix.

  • $\begingroup$ Accepting this as the answer because it told me what I needed to know. I computed the QR decomposition and checked the relative residual. For 32x32 size matrices I get 6E-15 residual, for 128x128 I get 7E-14, and for 512x512 I get 5.4E-13. Didn't have the patience to test bigger matrices. $\endgroup$ Jun 2 '20 at 13:22

Back-substitution is backward stable, and this implies that its error is of the same order as the error induced by truncation of the input data to machine precision, (machine precision)x(condition number of the problem).

So if back-substitution gives you an inaccurate answer, then no other algorithm can ever give you a better answer on the same problem (for generic data, apart from a worst-case factor O(n), typically even lower).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.