概要
よく知られているように，非可換の行列を同時対角化することはできません．しかしながら，「近似的」に同時対角化することは可能です．本稿では非可換行列を近似的に同時対角化する方法を紹介したいと思います．
方法
方法は文献の方法に倣います．この文献では，対角化したい行列群を $A^{(k)}, k=1,K$と表記するとして，$\sum_{k}\sum_{i<j}(U^H A^{(k)} U)_{ij}^2$という量が最小化される，すなわち作用した結果得られる行列の非対角項の和が最小化されるようなユニタリ変換$U$を求めることによって近似的に複数の行列を対角化します．
この目的のため，まずは以下のような回転行列を定義します．
R_{ij} =
\begin{pmatrix}
r_{ii} & r_{ij} \\
r_{ji} & r_{jj}
\end{pmatrix}
=
\begin{pmatrix}
c & s^\ast \\
-s & c^\ast
\end{pmatrix}
, c,s \in \mathbf{C}, |c|^2+|s|^2=1
この行列の行列要素は$ii, ij, ji, jj$以外の要素は単位行列と等しいです．この回転行列を作用させた関数$O$
$$
O=\sum_k \sum_{i<j}\left( R A^{(k)} R^\ast \right)_{ij}^2
$$
を最小化する$c,s$ペアを決めることができれば目的が達成できます．このため，次のような$3\times3$行列を定義します．
$$
G\equiv \Re\left(\sum_k h^\ast \left(A_k\right) h\left(A_k\right) \right), \\
h\left(A\right) \equiv \left[ A_{ii} - A_{jj}, A_{ij} + A_{ji}, i \left(A_{ji} - A_{ij} \right)\right ]
$$
目的の$c,s$ペアは$G$の最大固有値の固有ベクトル$x, y, z$を用いて以下のように求めることができます（その証明は文献を参照してください）
$$
c = \sqrt{\frac{x+r}{2r}}, s=\frac{y-iz}{\sqrt{2r\left(x+r\right)}},r=\sqrt{x^2+y^2+z^2}
$$
この手続きをすべての$ij$ペアで作用させて関数$O$を評価します．この処理を繰り返し行い，$O$が閾値の範囲で変化しなくなったら「最小化」できたとみなすことにします．もちろん制約の範囲での最小化なので，非対角項の和が絶対値として小さいとは限りません．回転行列による変換はユニタリ変換になっているため，何度作用しても結果得られる変換はユニタリ変換です．
実装例
ソースコード
Fortran90で非可換行列を同時対角化するプログラムを試作しました．少し長いですが，以下に全文掲載します．
program simultaneous_diagonalization_test
implicit none
integer :: ndim
real(kind=8) :: diag_factor
real(kind=8) :: nondiag_factor
integer :: nrand
real(kind=8) :: max_fluc
complex(kind=8), allocatable, dimension(:,:) :: testmatrix
complex(kind=8), allocatable, dimension(:,:,:) :: testmatrix_mult
real(kind=8), allocatable, dimension(:) :: eigenvalues,rwork
real(kind=8), allocatable, dimension(:,:) :: eigenvalues_mult
real(kind=8), allocatable, dimension(:,:) :: eigenvectors
complex(kind=8), allocatable, dimension(:) :: work
integer :: lwork,info,i,j,k,ierror
real(kind=8) :: rand
write(6,*) 'enter ndim'
read(5,*) ndim
write(6,*) 'enter diagonal and non-diagonal factor'
read(5,*) diag_factor,nondiag_factor
lwork = 2*ndim-1
allocate(testmatrix(ndim,ndim))
allocate(eigenvalues(ndim))
allocate(eigenvectors(ndim,ndim))
allocate(work(lwork))
allocate(rwork(3*ndim-2))
call get_test_matrix(ndim,diag_factor,nondiag_factor,testmatrix)
call zheev('V','U',ndim,testmatrix,ndim,eigenvalues,work,lwork,rwork,info)
write(6,*) 'eigenvalues from zheev'
do i=1,ndim
write(6,'(a,i8,f20.10)') 'eigenvalue ',i,eigenvalues(i)
enddo
allocate(testmatrix_mult(ndim,ndim,1))
call get_test_matrix(ndim,diag_factor,nondiag_factor,testmatrix)
testmatrix_mult(:,:,1) = testmatrix(:,:)
allocate(eigenvalues_mult(ndim,1))
call simultaneous_diagonalization(1,ndim,testmatrix_mult,eigenvalues_mult,eigenvectors,ierror)
write(6,*) 'eigenvalues from simultaneous_diagonalization'
do i=1,ndim
write(6,'(a,i8,f20.10)') 'eigenvalue ',i,eigenvalues_mult(i,1)
enddo
write(6,*) 'number of random matrices'
read(5,*) nrand
write(6,*) 'max. fluctuation'
read(5,*) max_fluc
deallocate(eigenvalues_mult)
deallocate(testmatrix_mult)
allocate(eigenvalues_mult(ndim,nrand))
allocate(testmatrix_mult(ndim,ndim,nrand))
do k=1,nrand
call get_test_matrix(ndim,diag_factor,nondiag_factor,testmatrix)
do i=1,ndim
do j=i,ndim
call random_number(rand)
testmatrix(i,j) = testmatrix(i,j)+(0.5d0-rand)*max_fluc*2
if(i/=j) testmatrix(j,i) = testmatrix(i,j)
enddo
enddo
testmatrix_mult(:,:,k) = testmatrix(:,:)
enddo
call simultaneous_diagonalization(nrand,ndim,testmatrix_mult,eigenvalues_mult,eigenvectors,ierror)
do k=1,nrand
write(6,*) 'eigenvalues for matrix ',k
do i=1,ndim
write(6,'(a,i8,f20.10)') 'eigenvalue ',i,eigenvalues_mult(i,k)
enddo
enddo
deallocate(eigenvalues)
deallocate(testmatrix)
deallocate(testmatrix_mult)
deallocate(eigenvalues_mult)
deallocate(work)
deallocate(rwork)
end program simultaneous_diagonalization_test
subroutine get_test_matrix(ndim,diag_factor,nondiag_factor,matrix)
implicit none
integer, intent(in) :: ndim
real(kind=8), intent(in) :: diag_factor,nondiag_factor
complex(kind=8), intent(out), dimension(ndim,ndim) :: matrix
integer :: i,j
do i=1,ndim
do j=i,ndim
if (i==j) then
matrix(i,j) = dcmplx(-diag_factor,0.d0)
else
matrix(i,j) = dcmplx(-nondiag_factor/(i-j)**2,0.d0)
matrix(j,i) = matrix(i,j)
endif
enddo
enddo
return
end subroutine get_test_matrix
function get_sum_offdiag(ndim, nmat, matrices) result(res)
implicit none
integer, intent(in) :: ndim,nmat
complex(kind=8), dimension(ndim,ndim,nmat) :: matrices
real(kind=8) :: res
integer :: i, j, k
res = 0.d0
do k=1, nmat
do i = 1,ndim-1
do j=i+1, ndim
res = res + real(dconjg(matrices(i,j,k))*matrices(i,j,k))
enddo
enddo
enddo
return
end function get_sum_offdiag
function get_sum_diag(ndim, nmat, matrices) result(res)
implicit none
integer, intent(in) :: ndim,nmat
complex(kind=8), dimension(ndim,ndim,nmat) :: matrices
real(kind=8) :: res
integer :: i, j, k
res = 0.d0
do k=1, nmat
do i = 1,ndim
res = res + real(dconjg(matrices(i,i,k))*matrices(i,i,k))
enddo
enddo
return
end function get_sum_diag
subroutine get_rotation_matrix_elements(ndim,nmat,matrices,i,j,rii,rjj,rij,rji)
implicit none
integer, intent(in) :: ndim,nmat
complex(kind=8), intent(in), dimension(ndim,ndim,nmat) :: matrices
integer, intent(in) :: i,j
complex(kind=8), intent(out) :: rii,rjj,rij,rji
real(kind=8), dimension(3,3) :: G
complex(kind=8), dimension(3) :: h
real(kind=8), dimension(3) :: eigenvalue_G
real(kind=8), dimension(3,3) :: hmat
complex(kind=8) :: aii,ajj,aij,aji
complex(kind=8) :: c,s
real(kind=8) :: x,y,z
complex(kind=8), parameter :: onei=(0.d0,1.d0)
integer, parameter :: lwork = 8
real(kind=8), dimension(lwork) :: work
integer :: k,inf,ii,jj
G = 0.d0
do k=1,nmat
aii=matrices(i,i,k)
ajj=matrices(j,j,k)
aij=matrices(i,j,k)
aji=matrices(j,i,k)
h(1) = aii-ajj
h(2) = aij+aji
h(3) = onei*(aji-aij)
do ii=1,3
do jj=1,3
hmat(ii,jj) = dconjg(h(ii))*h(jj)
enddo
enddo
G(:,:) = G(:,:)+real(hmat(:,:))
enddo
call dsyev('V','U',3,G,3,eigenvalue_G,work,lwork,inf)
if(G(1,3)<0) G(:,:) = -G(:,:)
x = G(1,3)
y = G(2,3)
z = G(3,3)
c = dcmplx(dsqrt(0.5d0+0.5d0*x),0.d0)
s = dcmplx(0.5d0*y,-0.5d0*z)/c
rii = c
rji = s
rij = -dconjg(s)
rjj = c
return
end subroutine get_rotation_matrix_elements
subroutine simple_sort(ndim,matrix,indices)
implicit none
integer, intent(in) :: ndim
real(kind=8), dimension(ndim), intent(inout) :: matrix
integer, dimension(ndim), intent(out) :: indices
real(kind=8) :: m
integer :: i, j, ii
do i=1,ndim
indices(i) = i
enddo
do i=1,ndim-1
do j=i+1,ndim
if(matrix(i)>matrix(j)) then
m = matrix(i)
matrix(i) = matrix(j)
matrix(j) = m
ii = indices(i)
indices(i) = indices(j)
indices(j) = ii
endif
enddo
enddo
return
end subroutine simple_sort
subroutine simultaneous_diagonalization(nmat,ndim,matrices,eigenvalues,eigenvectors,ierror)
implicit none
integer, intent(in) :: nmat, ndim
complex(kind=8), dimension(ndim,ndim,nmat), intent(inout) :: matrices
real(kind=8), dimension(ndim,nmat), intent(out) :: eigenvalues
real(kind=8), dimension(ndim,ndim), intent(out) :: eigenvectors
integer, intent(out) :: ierror
integer, parameter :: nsweep = 1000
real(kind=8), parameter :: eps = 1.d-8
complex(kind=8), parameter :: zero = (0.d0,0.d0)
complex(kind=8), parameter :: one = (1.d0,0.d0)
real(kind=8) :: sum_offdiag, sum_offdiag_prev
real(kind=8) :: sum_diag
complex(kind=8) :: rii,rjj,rij,rji,mi,mj
real(kind=8) :: get_sum_offdiag, get_sum_diag
complex(kind=8), allocatable, dimension(:) :: complex_i, complex_j
real(kind=8), allocatable, dimension(:) :: eigenvalue
integer, allocatable, dimension(:) :: indices
real(kind=8), allocatable, dimension(:,:) :: eigenvectors_buf
integer :: i, j, k, i0, isweep
sum_offdiag_prev = get_sum_offdiag(ndim, nmat, matrices)
eigenvectors = zero
do i=1,ndim
eigenvectors(i,i) = one
enddo
allocate(complex_i(ndim));allocate(complex_j(ndim))
do isweep = 1, nsweep
do i=1, ndim
do j=1, ndim
if(i==j) cycle
call get_rotation_matrix_elements(ndim,nmat,matrices,i,j,rii,rjj,rij,rji)
do k=1, nmat
do i0=1,ndim
mi = matrices(i,i0,k)
mj = matrices(j,i0,k)
matrices(i,i0,k) = mi*dconjg(rii)+mj*dconjg(rji)
matrices(j,i0,k) = mj*dconjg(rjj)+mi*dconjg(rij)
enddo
do i0=1,ndim
mi = matrices(i0,i,k)
mj = matrices(i0,j,k)
matrices(i0,i,k) = (rii)*mi+(rji)*mj
matrices(i0,j,k) = (rjj)*mj+(rij)*mi
enddo
enddo
do i0=1,ndim
complex_i(i0) = eigenvectors(i0,i)*rii+eigenvectors(i0,j)*rji
complex_j(i0) = eigenvectors(i0,i)*rij+eigenvectors(i0,j)*rjj
enddo
eigenvectors(:,i) = complex_i(:)
eigenvectors(:,j) = complex_j(:)
enddo
enddo
sum_offdiag = get_sum_offdiag(ndim, nmat, matrices)
sum_diag = get_sum_diag(ndim, nmat, matrices)
if(abs(sum_offdiag-sum_offdiag_prev)<eps) then
write(6,'(a,i6,a)') 'simultaneous diagonalization converged after ',isweep,' sweeps'
write(6,'(a,f20.15)') 'delta for the sum of the off-diagonal elements',abs(sum_offdiag-sum_offdiag_prev)
write(6,'(a,2f20.15)') 'abs. average value for the sum of the diagonal and offdiagonal elements ', &
abs(sum_diag)/dble(nmat)/dble(ndim),abs(sum_offdiag)/dble(nmat)/dble(ndim*(ndim-1)*0.5)
allocate(eigenvalue(ndim))
allocate(indices(ndim))
allocate(eigenvectors_buf(ndim,ndim))
do k=1,nmat
do i=1,ndim
eigenvalue(i) = real(matrices(i,i,k))
enddo
call simple_sort(ndim, eigenvalue, indices)
eigenvalues(:,k) = eigenvalue(:)
do i=1,ndim
eigenvectors(:,indices(i)) = eigenvectors_buf(:,i)
enddo
enddo
deallocate(complex_i)
deallocate(complex_j)
deallocate(eigenvalue)
deallocate(indices)
deallocate(eigenvectors_buf)
ierror = 0
return
endif
sum_offdiag_prev = sum_offdiag
enddo
write(6,'(a)') 'simultaneous diagonalization failed to converge'
ierror = 1
deallocate(complex_i)
deallocate(complex_j)
return
end subroutine simultaneous_diagonalization
- simultaneous_diagonalizationに同時対角化を行う処理が記述されています．
- get__rotation_matrix_elements で$r_{ii}, r_{ij}, r_{ji}, r_{jj}$の計算を行います．
- get_sum_diag, get_sum_offdiagでそれぞれ対角要素および非対角要素の和の計算を行います．
- get_test_matrixでテスト計算用の行列を作ります．テスト用の行列としては，はToeplitz行列を採用しました．対角要素は定数 $a$，非対角要素は$b/(i-j)^2$という具合に評価します．
- simple_sortで固有値を昇順にソートします．
- mainプログラム(サブルーチンsimultaneous_diagonalization_test)はまず単一の行列についてLAPACKとsimultaneous_diagonalizationで対角化します．ついで行列要素をランダムに変化させた複数の行列を同時対角化します．
コンパイルの仕方
コンパイルするにはLAPACKが必要です．たとえば以下の要領でコンパイルできるはずです．
$ gfortran simultaneous_diagonalization.F90 -llapack -lblas
使い方とテスト計算
実行するといろいろとパラメーターの指定を促されるので適宜入力します．
$ ./a.out
enter ndim # 行列の次元を入力する
10
enter diagonal and non-diagonal factor #対角及び非対角項のファクターを入力する
10.2 7.8
eigenvalues from zheev
eigenvalue 1 -30.7913801249
eigenvalue 2 -24.3381478761
eigenvalue 3 -18.6973305976
eigenvalue 4 -13.6783668636
eigenvalue 5 -9.3535576778
eigenvalue 6 -5.6854290655
eigenvalue 7 -2.6921957800
eigenvalue 8 -0.3619712059
eigenvalue 9 1.3003175438
eigenvalue 10 2.2980616475
simultaneous diagonalization converged after 5 sweeps
delta for the sum of the off-diagonal elements 0.000000000000000
abs. average value for the sum of the diagonal and offdiagonal elements 221.130642277836756 0.000000000000000
eigenvalues from simultaneous_diagonalization
eigenvalue 1 -30.7913801249
eigenvalue 2 -24.3381478761
eigenvalue 3 -18.6973305976
eigenvalue 4 -13.6783668636
eigenvalue 5 -9.3535576778
eigenvalue 6 -5.6854290655
eigenvalue 7 -2.6921957800
eigenvalue 8 -0.3619712059
eigenvalue 9 1.3003175438
eigenvalue 10 2.2980616475
number of random matrices # ランダムに値を追加する行列の数を指定します
4
max. fluctuation # ランダム変位の最大量を指定します
1.2
simultaneous diagonalization converged after 6 sweeps
delta for the sum of the off-diagonal elements 0.000000000330658
abs. average value for the sum of the diagonal and offdiagonal elements 223.916541066975071 0.340779305554451
eigenvalues for matrix 1
eigenvalue 1 -32.4113120479
eigenvalue 2 -25.0655263639
eigenvalue 3 -18.1536168634
eigenvalue 4 -12.3514343216
eigenvalue 5 -10.0010023358
eigenvalue 6 -4.7210201305
eigenvalue 7 -1.1943172339
eigenvalue 8 -0.3739783008
eigenvalue 9 0.4442298477
eigenvalue 10 1.7845332192
eigenvalues for matrix 2
eigenvalue 1 -30.1545092207
eigenvalue 2 -23.0262617720
eigenvalue 3 -20.1872100706
eigenvalue 4 -13.9406726704
eigenvalue 5 -9.4098994557
eigenvalue 6 -6.5666914451
eigenvalue 7 -1.6452553290
eigenvalue 8 -0.9816935756
eigenvalue 9 1.5071539917
eigenvalue 10 1.7755927175
eigenvalues for matrix 3
eigenvalue 1 -32.3851467052
eigenvalue 2 -22.9979774119
eigenvalue 3 -19.3975422400
eigenvalue 4 -13.0408574607
eigenvalue 5 -10.5791614500
eigenvalue 6 -5.1221915649
eigenvalue 7 -3.8254516664
eigenvalue 8 -0.8200685091
eigenvalue 9 2.5716528621
eigenvalue 10 3.4651432239
eigenvalues for matrix 4
eigenvalue 1 -30.6818354869
eigenvalue 2 -23.5782510883
eigenvalue 3 -19.3287057035
eigenvalue 4 -13.5552190879
eigenvalue 5 -9.0292689505
eigenvalue 6 -5.6978821278
eigenvalue 7 -3.6142706907
eigenvalue 8 1.0435680522
eigenvalue 9 1.6057939078
eigenvalue 10 2.2128183463
入力パラメーターと結果をそのまま掲載してみました．まず，次元が10, 対角要素のファクターが10.2, 非対角項のファクターが7.8という設定でテスト用行列を作りました．この行列の固有値の計算結果はLAPACKとsimultaneous_diagonalizationとで完全に一致しました．非対角項の和もゼロとなりよい具合です．
つぎに，もとの行列に対し各行列要素が最大1.2だけランダムに変位した行列を4つ用意し，同時対角化を実行しました．今度は非対角項はゼロにはならず，0.340779305554451という値になりました．これは対角項の和の223.916541066975071と比較すると小さく，目論見通り近似的な同時対角化に成功していそうです．
計算速度
ここで採用した手法は，$RA^{(k)}R^*$の行列－行列積が$O(N)$の演算量で，この演算をすべての$i,j$ペアで行う必要があるので$O(N^3)$に比例する演算量となります（実対称/エルミート行列の場合はそれ相応に演算量を減らせますが，試作プログラムではそのことを考慮していません）．具体的に行列のサイズを変えて計算時間を検証してみました．あわせてLAPACK(最適化されていないNetlib版)の計算時間も示します．単位は秒です．
| 行列サイズ | LAPACKの処理時間 | simultaneous_diagonalizationの処理時間 |
|---|---|---|
| 10 | 1.3e-4 | 1.2e-3 |
| 50 | 2.6e-3 | 8.6e-2 |
| 100 | 1.7e-2 | 0.5 |
| 500 | 0.94 | 53.2 |
| 1000 | 5.7 | 435.9 |
行列サイズと計算時間の関係から，おおよそ$O(N^3)$に比例する処理時間がかかっています．$f(x) = ax^b$のような関数にフィットすると，いずれの場合も$b=2.9$程度となりました．LAPACKの方がはるかに高速です．おそらくチューニングなどで埋められる差ではなく，単一の行列を対角化するのならば伝統的な方法を使った方がよさそうです．
参考文献
Jean- François Cardoso and Antoine Souloumiac, “Jacobi Angles for Simultaneous Diagonalization”, Journal of Mathematical Analysis and Applications 17 161-164 (1996).