fortranspatial-interpolationqr-decomposition

QR decomposition with increasingly large errors


I received a code to work on and it has a QR decomposition subroutine as shown below. When checking this QR decomposition using the error calculation I added at the end of it, I found the error increases to very high magnitudes after hundreds of simulation steps.

The QR decomposition is used to solve the coefficients used in triquadratic interpolation, adapted from trilinear interpolation as found on Wikipedia: https://en.wikipedia.org/wiki/Trilinear_interpolation. In the code, we generate a 5x5x5 grid of values, and uses triquadratic interpolation to derive a spatial function form of the values to obtain f(x,y,z) = a_0 + a_1x + a_2y + a_3*z + ...

The function f(x,y,z) can be written as Ax = b for the pre-existing values on each grid point, such that b is the 125x1 vector of existing values, A is a 125x20 matrix calculated using the coordinates, and x is the 20x1 coefficient vector composed of (a_0, a_1, a_2, ...).

Then, x is calculated using the QR decomposition algorithm shown below. When I looked up QR decomposition, I found some discussions here that mentioned certain approaches can become increasingly unstable, but approaches such as householder reflection should be very stable.

My first question: I am having trouble identifying the type of approach used in my code, it doesn't normalize the variables and I can't match the process to other implementations found online. Can anyone tell me if this implementation is already a householder reflection?

My second question: If the implementation is correct and should be stable, can anyone give me some ideas to check why the solution would become increasingly unstable? Or if it is not householder or simply a bad implementation, would an example householder implementation that I found online help improve the accuracy?

  ! QR decomposition solving for vector B in AX = B
  subroutine qr_eqsystem(A, B, X, M, N, error)
  
  implicit none
  
  integer, intent(in) :: m, n
  real, intent(in) :: A(m, n), B(m)
  
  real, intent(out) :: x(n)
  real, intent(out) :: error
  
  !locals
  real :: q(m, n), r(n, n), aux
  integer :: i, j, k
  
  do i = 1,n
    do k = 1,m
      q(k, i) = A(k, i)
    enddo
  enddo
  
  do i = 1,n
    do j = 1,n
      r(j, i) = 0.0
    enddo
  enddo
  
  do i = 1,n
    do j = 1,i-1
      r(j, i) = 0.0
      do k = 1,m
        r(j, i) = r(j, i) + q(k, j)*A(k, i)
      enddo
      
      do k = 1,m
        q(k, i) = q(k, i) - r(j, i)*q(k,j)
      enddo
    enddo
    
    r(i, i) = 0.0
    do k = 1,m
      r(i, i) = r(i, i) + q(k, i)*q(k, i)
    enddo
    r(i, i) = sqrt(r(i, i))
    
    if(r(i,i).ne.0.0)then
      do k = 1,m
        q(k, i) = q(k, i) / r(i, i)
      enddo
    endif
    
  enddo
  
  do i = 1,n
    x(i) = 0.0
    do j = 1,m
      x(i) = x(i) + q(j, i)*B(j)
    enddo
  enddo
  
  if(r(n, n).ne.0.0)then
    x(n) = x(n) / r(n,n)
  endif
  
  do i = n-1,1,-1
    do j = i+1,n
      x(i) = x(i) - r(i, j) * x(j)
    enddo
    
    if(r(i, i).ne.0.0)then
       x(i) = x(i) / r(i, i)
    endif
  enddo
  
  error = 0
  
  ! calculate AX and compare with B as an error, report total raw error
  do i = 1,m
    aux = 0.0
    do j = 1,n
      aux = aux + A(i,j)*X(j)
    enddo
    error = error + abs(aux - B(i))
  enddo
  
  endsubroutine qr_eqsystem

Solution

  • Your QR decomposition implementation looks normal (far from ideal, but normal). You can compare the absolute values of the diagonal R with the LAPACK implementation or with Python numpy.linalg.qr().

    However, numerical stability of the solution only means that

    |x - x'|/|x| <= f(m, n)*cond(A)*eps

    where, f(m, n) is a function of a small degree from m and n, cond(A) is a condition number approximately equal to the ratio of the maximum and minimum absolute values diagonal of the matrix R (A ≈ QR).

    However, the magnitude of the residual of the solution of the least squares system Ax = B may still be significant for two reasons:

    1. The system is incompatible;
    2. The number of conditionality (cond(A)) is too high.

    I think it's worth looking for computational errors in calculating the matrix A and the right parts of B.