QR-разложение прошло успешно, но имеет большое количество ошибок
У меня есть две реализации QR-декомпозиции в коде на Фортране. Одна из них — это пользовательская функция разложения QR:
! QR decomposition solving for vector B in AX = B
subroutine qr_eqsystem(A, B, X, M, N, q, r, error)
implicit none
!CALL qr_eqsystem(A, C(:,i,j), B, GNDINT, NEQ, IERR)
integer, intent(in) :: m, n
real, intent(in) :: A(m, n), B(m)
real, intent(out) :: x(n),q(m, n),r(n, n)
real, intent(out) :: error
!locals
real :: 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
error = sum(abs(matmul(A,X)-B))
endsubroutine qr_eqsystem
Второй использует встроенный LAPACK Intel Fortran с dgeqrf, dormqr и dtrtrs:
q = A
call dgeqrf(totPoints,ncol,q,totPoints,tau,work,1280,info)
Call dormqr('L','T',totPoints,1,ncol,q,totPoints,tau,C,
# totPoints,work,1280,info)
x = C
Call dtrtrs('U','N','N',ncol,1,q,totPoints,x,totPoints,info)
error = sum(abs(matmul(A,x(1:ncol))-C))
Приведенные выше коды позволяют найти x в системе уравнений Ax = B, где размер A равен 125x21, B 125x1 и x 21x1. Во время моделирования матрица A является относительно постоянной, так как она построена на основе пространственных координат, но вектор B развивается со временем, при этом величина элементов в B и разница между самым большим и наименьшим элементами со временем увеличиваются. Оба кода реализованы в коде Fortran для отладки, где пользовательская реализация является основным решателем, и когда ошибка превышает высокое значение, например 10,0, код вызывает подпрограмму Intel для проверки.
Значения error в обоих кодах возвращают очень низкие значения (<1e-5) в начале моделирования, но error увеличиваются до высоких значений по мере продолжения моделирования и, похоже, не имеют насыщенного значения, и я этого не делаю. понять, почему это происходит. Я проверил разложенные матрицы Q и R: Q^T*Q возвращает единичную матрицу, а Q*R возвращает A с разумной точностью. Насколько я понимаю подход к декомпозиции QR, поскольку A успешно разлагается, x всегда должен быть разрешимым, поскольку он включает в себя простое вычисление Q^T*B, а обратная замена Rx = Q^T*B для решения x всегда должна быть возможна, поскольку все значения действительны и мои переменные достаточно велики (8 байт с реальным размером флага: 64), чтобы не вызывать значительных ошибок с плавающей запятой.
Может ли кто-нибудь дать мне больше информации о том, почему QR-разложение Ax = B может дать большую ошибку при оценке A*xp - B, где xp является решенным x из процесса?
Кроме того, я проверил MATLAB x = A\B. Интересно, что решения Intel Fortran dgeqrf и MATLAB дают более высокие ошибки, чем первая пользовательская реализация, причем процесс MATLAB дает наибольшую ошибку, несмотря на то, что разложенные матрицы Q и R из всех трех подходов практически идентичны друг другу.
Код для матрицы A. Это написано в фиксированной форме Фортрана с расширением .for. GridCoords просто генерируется путем сопоставления сетки точек 5x5x5 в диапазоне [-1 1] с диапазоном, соответствующим объекту размером 1e-5, например [-2.5e-5 2.5e-5] для объекта диаметром 50 микрон.
subroutine defineA(ix,iy,iz,A,ipt)
implicit none
integer, intent(in) :: ix,iy,iz,ipt
real, intent(out) :: A(totPoints,ncol)
if (iInterp.ge.1) then
A(ipt,1) = 1
A(ipt,2) = gridCoords(1,ix)
A(ipt,3) = gridCoords(2,iy)
A(ipt,4) = gridCoords(3,iz)
A(ipt,5) = gridCoords(1,ix)*gridCoords(2,iy)
A(ipt,6) = gridCoords(1,ix)*gridCoords(3,iz)
A(ipt,7) = gridCoords(2,iy)*gridCoords(3,iz)
A(ipt,8) = gridCoords(1,ix)*gridCoords(2,iy)*gridCoords(3,iz)
if (iInterp.ge.2) then
A(ipt,9) = gridCoords(1,ix)*gridCoords(1,ix)
A(ipt,10) = gridCoords(2,iy)*gridCoords(2,iy)
A(ipt,11) = gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,12) = gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(2,iy)
A(ipt,13) = gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(3,iz)
A(ipt,14) = gridCoords(1,ix)*gridCoords(2,iy)*gridCoords(2,iy)
A(ipt,15) = gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(3,iz)
A(ipt,16) = gridCoords(1,ix)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,17) = gridCoords(2,iy)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,18) = gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)
A(ipt,19) = gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,20) = gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,21) = gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(3,iz)
A(ipt,22) = gridCoords(1,ix)*gridCoords(2,iy)
# *gridCoords(2,iy)*gridCoords(3,iz)
A(ipt,23) = gridCoords(1,ix)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,24) = gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)
A(ipt,25) = gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(3,iz)
# *gridCoords(3,iz)
A(ipt,26) = gridCoords(1,ix)*gridCoords(2,iy)
# *gridCoords(2,iy)*gridCoords(3,iz)
# *gridCoords(3,iz)
A(ipt,27) = gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)
if (iInterp.ge.3) then
A(ipt,28)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
A(ipt,29)=gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
A(ipt,30)=gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,31)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)
A(ipt,32)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(3,iz)
A(ipt,33)=gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
A(ipt,34)=gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)
A(ipt,35)=gridCoords(1,ix)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,36)=gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,37)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)
# *gridCoords(3,iz)
A(ipt,38)=gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)
A(ipt,39)=gridCoords(1,ix)
# *gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,40)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)
A(ipt,41)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,42)=gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
A(ipt,43)=gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,44)=gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,45)=gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,46)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)
A(ipt,47)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,48)=gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,49)=gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)
A(ipt,50)=gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,51)=gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,52)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,53)=gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,54)=gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,55)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
A(ipt,56)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,57)=gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,58)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)
A(ipt,59)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,60)=gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,61)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,62)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,63)=gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
A(ipt,64)=gridCoords(1,ix)*gridCoords(1,ix)*gridCoords(1,ix)
# *gridCoords(2,iy)*gridCoords(2,iy)*gridCoords(2,iy)
# *gridCoords(3,iz)*gridCoords(3,iz)*gridCoords(3,iz)
endif
endif
endif
endsubroutine defineA
@MartinBrown Я удалил свой предыдущий комментарий. Симуляция достигла той же точки, что и раньше, когда я установил условие ошибки > 10, и на том же шаге она все равно достигла ошибки > 10. Но на этот раз ошибка процесса dgeqrf имеет величину 3000, в то время как ошибка пользовательской реализации, использующей реальные 16, по-прежнему равна только 16. Теперь я прогоню всю симуляцию, чтобы увидеть ошибку в конце, и посмотрю еще несколько способов устранения ошибки с плавающей запятой. Возможно, я смогу масштабировать матрицу A, чтобы она не использовала такие маленькие значения.
Я думаю, стоит посмотреть, как вычисляются матрица A и вектор B и величины входящих в них переменных. Я подозреваю, что проблема как-то затруднительна по мере приближения к вашему решению. Вам необходимо улучшить его состояние.
Возможно, стоит опубликовать код, создающий матрицу A и вектор B. Вы производите выборку своей функции с равными шагами по x/y/z или делаете что-то более сложное?
@MartinBrown Матрица A проста, добавил я в пост. Вектор B гораздо сложнее. Он выбирается из распределения, рассчитанного на основе второго моментного разброса, который определяет стандартное отклонение, затем присваивается каждой из точек сетки 5x5x5, а 125 значений располагаются по столбцам в векторе B, соответствующем матрице A.
Я могу привести примеры вектора B, используемого в моем коде, но коды, которые его создали, более сложны, и я не думаю, что мне следует публиковать их в Интернете в любом случае.
Я думаю, что, по крайней мере, часть проблемы заключается в том, что ваше масштабирование переменных x,y,z~10^-5 вместе с задействованными высокими степенями делает матрицу A чрезвычайно близкой к единственному числу и подверженной большим ошибкам округления. Другой вопрос, который у меня возникает, заключается в том, как 5 известных точек данных распределены по [-1,1] -1 -0,5 0 0,5 1 или чему-то вроде нулей полинома Чебышева T5 (x). Поскольку это зависит только от координат x, y, z, почему A вообще меняется? (я ожидал, что это также будет зависеть от ваших известных данных). Также нормально ли это работает для iInterp = 1 & 2?
Ответы 1
У меня есть две реализации QR-разложения... Приведенные выше коды определяют x в системе уравнений Ax = B, где размер A равен 125x21, B равен 125x1, а x равен 21x1...
Обратите внимание, что количество уравнений в вашей системе превышает количество переменных, поэтому точного решения не существует. И решением вашей Ax = B-системы, результирующим решением Rx = Q^T*B-системы, будет такой x, для которого достигается минимальный функционал (B - Ax)*(B - Ax). Другими словами, это решение задачи наименьших квадратов.
Значения ошибки в обоих кодах возвращают очень низкие значения (<1e-5) в начале моделирования, но ошибка увеличивается до высоких величин по мере продолжения моделирования и, похоже, не имеет насыщенного значения, и я этого не делаю. понять, почему это происходит.
В предыдущем вопросе вы написали: «QR-разложение используется для решения коэффициентов, используемых в квадратичной интерполяции, адаптированной из трилинейной...». Я не очень понимаю, откуда у вас функция с 21 коэффициентом: f(x,y,z) = a_0 + a_1x + a_2y + a_3*z+ ....
Но вполне вероятно, что в процессе моделирования целевая функция стремится к виду, который не может быть хорошо представлен функцией f(x,y,z) с 21 коэффициентом. Таким образом, решение вашей системы начинает иметь большой остаток.
Кроме того, я проверил x = A\B MATLAB. Интересно, что и решения Intel Fortran dgeqrf, и решения MATLAB дают более высокие ошибки, чем первая пользовательская реализация, при этом процесс MATLAB дает наибольшую ошибку, несмотря на то, что разложенные матрицы Q и R из всех трех подходов практически идентичны друг другу.
Поскольку вы сравниваете error = sum(abs(matmul(A,X)-B)), но qr_eqsystem(), Intel Fortran dgeqrf и MATLAB, минимизируйте error2 = DOT_PRODUCT(matmul(A,X)-B, matmul(A,X)-B), то это сравнение ничего не говорит о качестве этих трёх решений. Вполне возможно, что более точное решение задачи наименьших квадратов (в норме L2) будет иметь худшую невязку в норме L1 (сумму абсолютных значений отклонений).
Если нет физической причины использовать норму Л1, лучше использовать норму Л2.
Во время моделирования матрица A является относительно постоянной, поскольку она строится из пространственных координат...
Если A — константа, то cond(A) = cond(R) ∼ max(abs(diagonal(R)))/min(abs(diagonal(R))) это тоже константа. Соответственно, если проблема решена хорошо в начале моделирования, то причина, скорее всего, не в точности матрицы, ее декомпозиции или решении.
Резюме: поскольку вы пишете «относительно постоянный», все равно стоит посмотреть на изменение cond(A) во время моделирования.
Если cond(A) > 10¹⁰ в начале или во время моделирования, то может возникнуть проблема с точностью решения и/или, возможно, проблема с координатами функции f(x,y,z);
Если cond(A) всегда < 10⁶, то проблема определенно в том, что форма вашего f(x,y,z) не соответствует модели. Тогда его нужно заменить на что-то другое, может на многомерные B-сплайны, может на что-то более подходящее к физике задачи.
Это очень похоже на предыдущий пост. Я подозреваю, что именно моделирование для создания подходящей матрицы
Aвне решателя является основной причиной ваших проблем. МатрицаA, должно быть, балансирует на грани нестабильности или что-то в этом роде. Заставить решатель матрицы работать в арифметике вещественного*10 или вещественного*16 может быть достаточно, чтобы взять его под контроль.