非可換行列を同時対角化する方法

概要

よく知られているように，非可換の行列を同時対角化することはできません．しかしながら，「近似的」に同時対角化することは可能です．本稿では非可換行列を近似的に同時対角化する方法を紹介したいと思います．

方法

方法は文献の方法に倣います．この文献では，対角化したい行列群を $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).