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
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:
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.