Алгоритм извлечения квадратного корня из матрицы

ulanna

Конкретно: из симметричной положительной матрицы с единицами на диагонали. Очень нужен этот алгоритм, чтобы потом на его основе написать макрос в Экселе. Очень, очень сильно нужен. Требований по скорости нет, точность - ну хоть какая-нибудь чтоб была. Потребуется для работы с матрицами размера до 20х20. Если есть у кого-нибудь указанный алгоритм в каком-либо виде - пожалуйста, поделитесь! В долгу не останусь!

Nitochka

Если бы не для экселя - то дал бы..
А так нуна искать собственные значения и векторы, извлекать квадратный корень из диагональной матрицы собственных значений, и делать обратное преобразование
A=O^{+}DO -> D^{1/2} -> A^{1/2}=O^{+}D^{1/2}O
Ну если очень хоца - могу дать фортрановские исходники по диагонализации матрицы. Тока вряд ли тебе это поможет

ulanna

Если бы не для экселя - то дал бы..
Да неважно. Я переведу. Можно в Си, в Паскале.
А так нуна искать собственные значения и векторы, извлекать квадратный корень из диагональной матрицы собственных значений, и делать обратное преобразование

Ну это-то всё понятно, но только я не знаю (точнее, не помню) как это всё делается. Причём, насколько я помню, осуществляется это довольно-таки непросто.
Так что, поможешь?

griz_a

Тебе нужен алгоритм разложения симметричной матрицы в произведение C^{-1}DC?

Nitochka

Ну держи:
Итак, tred2 сводит исходную симметричную (!) матрицу к трехдиагональной, tqli довершает процесс. pythag - вспомогательная хрень.
Заточены под работу последовательно. Да, оригинал tqli я куда-то проимел, эта слегка модифицированная версия под real*8 арифметику. Можно всегда найти в инете по ключевому слову Numerical Recipes
      
SUBROUTINE tqli(d,e,n,np,z)
С d - диагональ исходной матрицы, e - наддиагональ, n - размерность оной, np - размер массива под эту матрицу, z - матрица собственных векторов
C тонкость: e(1)=0, заполнение от e(2) до e(n)
INTEGER n,np
real*8 d(npe(npz(np,np)
CU USES pythag
INTEGER i,iter,k,l,m
real*8 b,c,dd,f,g,p,r,s,pythag
do 11 i=2,n
e(i-1)=e(i)
11 continue
e(n)=0.0d0
do 15 l=1,n
iter=0
1 do 12 m=l,n-1
dd=dabs(d(m+dabs(d(m+1
if (dabs(e(m+dd.eq.dd) goto 2
12 continue
m=n
2 if(m.ne.l)then
if(iter.eq.30)pause 'too many iterations in tqli'
C if(iter.eq.30) write(6,'(''too many iterations in tqli'')')

iter=iter+1
C*
C write(6,'(2x,''iter='',i4)') iter
C pause
C*
g=(d(l+1)-d(l/(2.0d0*e(l
r=pythag(g,1.0d0)
g=d(m)-d(l)+e(l)/(g+sign(r,g
s=1.0d0
c=1.0d0
p=0.0d0
do 14 i=m-1,l,-1
f=s*e(i)
b=c*e(i)
r=pythag(f,g)
e(i+1)=r
if(r.eq.0.0d0)then
d(i+1)=d(i+1)-p
e(m)=0.0d0
goto 1
endif
s=f/r
c=g/r
g=d(i+1)-p
r=(d(i)-g)*s+2.0d0*c*b
p=s*r
d(i+1)=g+p
g=c*r-b
C Omit lines from here ...
do 13 k=1,n
f=z(k,i+1)
z(k,i+1)=s*z(k,i)+c*f
z(k,i)=c*z(k,i)-s*f
13 continue
C ... to here when finding only eigenvalues.
14 continue
d(l)=d(l)-p
e(l)=g
e(m)=0.0d0
goto 1
endif
15 continue
return
END

FUNCTION pythag(a,b)
real*8 a,b,pythag
real*8 absa,absb
absa=dabs(a)
absb=dabs(b)
if(absa.gt.absb)then
pythag=absa*dsqrt(1.0d0+(absb/absa)**2)
else
if(absb.eq.0.0d0)then
pythag=0.0d0
else
pythag=absb*dsqrt(1.0d0+(absa/absb)**2)
endif
endif
return
END

SUBROUTINE tred2(a,n,np,d,e)
INTEGER n,np
REAL*8 a(np,npd(npe(np)
!Householder reduction of a real, symmetric, n by n matrix a,
!stored in an np by np physical array. On output, a is replaced
!by the orthogonal matrix Q e?ecting the transformation. d
!returns the diagonal elements of the tridiagonal matrix, and
!e the o?-diagonal elements, with e(1)=0. Several statements,
!as noted in comments, can be omitted if only eigenvalues are
!to be found, in which case a contains no useful information
!on output. Otherwise they are to be included.
INTEGER i,j,k,l
REAL*8 f,g,h,hh,scale
do i=n,2,-1
l=i-1
h=0.d0
scale=0.d0
if(l.gt.1)then
do k=1,l
scale=scale+abs(a(i,k
enddo
if(scale.eq.0.)then !Skip transformation.
e(i)=a(i,l)
else
do k=1,l
a(i,k)=a(i,k)/scale !Use scaled aТs for transformation.
h=h+a(i,k)**2 !Form ? in h.
enddo
f=a(i,l)
g=-sign(dsqrt(hf)
e(i)=scale*g
h=h-f*g !Now h is equation (11.2.4).
a(i,l)=f-g !Store u in the ith row of a.
f=0.d0
do j=1,l
C Omit following line if finding only eigenvalues
a(j,i)=a(i,j)/h !Store u/H in ith column of a.
g=0.d0 !Form an element of A ╖ u in g.
do k=1,j
g=g+a(j,k)*a(i,k)
enddo
do k=j+1,l
g=g+a(k,j)*a(i,k)
enddo
e(j)=g/h !Form element of p in temporarily unused
f=f+e(j)*a(i,j) !element of e.
enddo
hh=f/(h+h) !form K, equation (11.2.11).
do j=1,l !Form q and store in e overwriting p.
f=a(i,j)
g=e(j)-hh*f
e(j)=g
do k=1,j !Reduce a, equation (11.2.13).
a(j,k)=a(j,k)-f*e(k)-g*a(i,k)
enddo
enddo
endif
else
e(i)=a(i,l)
endif
d(i)=h
enddo
C Omit following line if finding only eigenvalues.
d(1)=0.d0
e(1)=0.d0
do i=1,n !Begin accumulation of transformation matrices.
C Delete lines from here ...
l=i-1
if(d(i).ne.0.d0)then !This block skipped when i=1.
do j=1,l
g=0.d0
do k=1,l !Use u and u/H stored in a to form P ╖ Q.
g=g+a(i,k)*a(k,j)
enddo
do k=1,l
a(k,j)=a(k,j)-g*a(k,i)
enddo
enddo
endif
C ... to here when finding only eigenvalues.
d(i)=a(i,i) !This statement remains.
C Also delete lines from here ...
a(i,i)=1.d0 !Reset row and column of a to identity matrix for
!next iteration.
do j=1,l
a(i,j)=0.d0
a(j,i)=0.d0
enddo
C ... to here when finding only eigenvalues.
enddo
return
END

Ну а уж взятие квадратного корня и обратное преобразование - сам напишешь

ulanna

Тебе нужен алгоритм разложения симметричной матрицы в произведение C^{-1}DC?
Ага!

ulanna

Офигеть...

Nitochka

Чо, впечатлило?
Стандартные процедуры из Numerical Recipes - реализация алгоритма QR/QL

MammonoK

Denman-Beavers square root iteration
Y_0 = A, Z_0 = I
Y_{k+1} = (Y_k + Z_k^{-1}) / 2
Z_{k+1} = (Z_k + Y_k^{-1}) / 2,
тогда
\lim Y_k = A^{1/2}
А вот для symmetric positive definite matrix:
Given a symmetric positive definite matrix A \in R^{n \times n}
this algorithm computes X = A^{1/2}
1. A = R^{T} R (Cholesky factorization)
2. Compute U = X_{\infty} from X_{k+1} = 1/2 (X_k + X_k^{-T}, k = 1, 2, ...
with X_0 = R.
3. X = U^{T}R.
В первом алгоритме фактически надо только реализовать нахождение обратной матрицы.

Nitochka

Второй алгоритм - практически один-в-один QR/QL

ulanna

Denman-Beavers square root iteration
Y_0 = A, Z_0 = I
Y_{k+1} = (Y_k + Z_k^{-1}) / 2
Z_{k+1} = (Z_k + Y_k^{-1}) / 2,
тогда
\lim Y_k = A^{1/2}
Во, вот это я понимаю. Обратную матрицу найти смогу
Оставить комментарий
Имя или ник:
Комментарий: