Ad と ad - fortran66のブログ

行列の指数関数は、べき展開した指数関数を定義として行列積であらわされます。この場合も収束半径は∞になるので、（オーバーフローしなければ）力づくで計算すればそれなりに結果が出てきます。

Ad は、演算子の同値変換を表わし行列であらわした場合、Ad(X)Y = XYX^-1 で定義されます。この変換の無限小変換は ad(x)Y = [x, Y] = xY - Yx と交換子積になります。ここでx は X の無限小変換で、X = exp(x) の関係があります。ad(x) は演算子で、ad(x)^2Y=ad(x)(ad(x)Y)=ad(x)[x, Y]=[x,[x, Y]] のように作用します。

また無限小変換を積分することで Ad(X)Y=exp(ad(x)Y) の関係が導かれます。　

ここでは X=exp(-itx) として、上記の関係式を数値的に確かめてみます。この場合、左辺は Ad(X)Y = exp(-itx) Y exp(itx) となり、exp(itx) と exp(-itx) の二つの行列の指数関数を求めて、その後 Y と行列積を計算します。一方、右辺は、exp(-itad(x)Y) を計算しなければなりませんが、ad(x) が演算子であるため、指数関数のべきのところで交換子積をとる専用の関数を書いて計算します。

右辺と左辺の計算値を書き出して、その一致をみます。なお、ここで行列は複素数にとって、三次元空間の回転行列を例に X軸周りと Y軸周りの回転を演算子としてとります。回転角を変えてみて、どの角度でも行列が一致することを確かめます。

exp 関数はオーバーロード定義しました。（Fortran のデフォルトでは、要素毎の exp を取る。）

ソースプログラム

    module m_matexp
      implicit none
      interface exp 
        procedure :: cmatexp
      end interface
    contains
      function commutator(x, y) result(res) ! [x, y]
        complex, intent(in) :: x(:, :),  y(:, :)
        complex :: res(size(x, 1), size(x, 2))
        res = matmul(x, y) - matmul(y, x)
      end function commutator
    
      function cmatexp(x) result(res)
        complex, intent(in) :: x(:, :)
        complex :: res(size(x, 1), size(x, 2)), tmp(size(x, 1), size(x, 2))
        integer :: i
        tmp = (0.0, 0.0)
        res = (0.0, 0.0)
        forall(i = 1:size(res, 1)) tmp(i, i) = (1.0, 0.0)
        do i = 1, 1000
          res = res + tmp   
          tmp = matmul(tmp, x) / i 
          if (maxval(abs(tmp)) < epsilon(0.0)) exit 
        end do
      end function cmatexp
      
      function Adj(x, Y) result(res) ! exp(x)Yexp(-x)
        complex, intent(in) :: x(:, :),  Y(:, :)
        complex :: res(size(x, 1), size(x, 2))
        res = matmul(matmul(exp(x), Y), exp(-x))
      end function Adj
      
      function ad_matexp(x, y) result(res)
        complex, intent(in) :: x(:, :), y(:, :)
        complex :: res(size(x, 1), size(x, 2)), tmp(size(x, 1), size(x, 2))
        integer :: i
        tmp = y
        res = (0.0, 0.0)
        do i = 1, 1000
          res = res + tmp   
          if (maxval(abs(tmp)) < epsilon(0.0)) exit 
          tmp = commutator(x, tmp) / i 
        end do
      end function ad_matexp
    end module m_matexp  
    
    program adjoint
      use m_matexp
      implicit none
      real, parameter :: pi = 4 * atan(1.0)
      integer, parameter :: n = 3
      complex :: lx(n, n), ly(n,n), lz(n, n)
      complex :: x(n, n), y(n,n), cy(n, n), cz(n, n)
      real :: theta
      integer :: i
      lx = transpose(reshape([[(0.0, 0.0), (1.0, 0.0), (0.0, 0.0)], &
                              [(1.0, 0.0), (0.0, 0.0), (1.0, 0.0)], &
                              [(0.0, 0.0), (1.0, 0.0), (0.0, 0.0)]], [n, n]) / sqrt(2.0))       
      ly = transpose(reshape([[(0.0, 0.0), (0.0,-1.0), (0.0, 0.0)], &
                              [(0.0, 1.0), (0.0, 0.0), (0.0,-1.0)], &
                              [(0.0, 0.0), (0.0, 1.0), (0.0, 0.0)]], [n, n]) / sqrt(2.0))       
      lz = transpose(reshape([[(1.0, 0.0), (0.0, 0.0), (0.0, 0.0)], &
                              [(0.0, 0.0), (0.0, 0.0), (0.0, 0.0)], &
                              [(0.0, 0.0), (0.0, 0.0), (-1.0, 0.0)]], [n, n]) )       
      
      do i = 1, 5 
        theta = pi / 4 * (i - 1)  
        x  = -theta * (0.0, 1.0) * lx
        Y = exp(-theta * (0.0, 1.0) * ly)

        cy = Adj(x, Y)
        print *, 't =', theta
        print *, 'Adjoint: Ad(itx)Y = exp(-itx) Y exp(itx)'
        print '(3("(", 2f8.5, ")"))', cy
      
        cz = ad_matexp(x, Y)
        print *, 'adjoint: exp(ad(itx))Y'
        print '(3("(", 2f8.5, ")"))', cz
        print *
      end do
    end program adjoint

実行結果

 t =  0.0000000E+00
 Adjoint: Ad(itx)Y = exp(-itx) Y exp(itx)
( 1.00000 0.00000)( 0.00000 0.00000)( 0.00000 0.00000)
( 0.00000 0.00000)( 1.00000 0.00000)( 0.00000 0.00000)
( 0.00000 0.00000)( 0.00000 0.00000)( 1.00000 0.00000)
 adjoint: exp(ad(itx))Y
( 1.00000 0.00000)( 0.00000 0.00000)( 0.00000 0.00000)
( 0.00000 0.00000)( 1.00000 0.00000)( 0.00000 0.00000)
( 0.00000 0.00000)( 0.00000 0.00000)( 1.00000 0.00000)

 t =  0.7853982
 Adjoint: Ad(itx)Y = exp(-itx) Y exp(itx)
( 0.78033-0.50000)( 0.35355-0.10355)( 0.07322 0.00000)
(-0.35355 0.10355)( 0.85355 0.00000)( 0.35355 0.10355)
( 0.07322 0.00000)(-0.35355-0.10355)( 0.78033 0.50000)
 adjoint: exp(ad(itx))Y
( 0.78033-0.50000)( 0.35355-0.10355)( 0.07322 0.00000)
(-0.35355 0.10355)( 0.85355 0.00000)( 0.35355 0.10355)
( 0.07322 0.00000)(-0.35355-0.10355)( 0.78033 0.50000)

 t =   1.570796
 Adjoint: Ad(itx)Y = exp(-itx) Y exp(itx)
(-0.00000-1.00000)(-0.00000 0.00000)( 0.00000-0.00000)
( 0.00000-0.00000)( 1.00000 0.00000)(-0.00000-0.00000)
( 0.00000 0.00000)( 0.00000 0.00000)(-0.00000 1.00000)
 adjoint: exp(ad(itx))Y
(-0.00000-1.00000)( 0.00000 0.00000)(-0.00000 0.00000)
(-0.00000-0.00000)( 1.00000 0.00000)( 0.00000-0.00000)
(-0.00000 0.00000)(-0.00000 0.00000)(-0.00000 1.00000)

 t =   2.356194
 Adjoint: Ad(itx)Y = exp(-itx) Y exp(itx)
(-0.28033-0.50000)(-0.35355 0.60355)( 0.42678 0.00000)
( 0.35355-0.60355)( 0.14645 0.00000)(-0.35355-0.60355)
( 0.42678 0.00000)( 0.35355 0.60355)(-0.28033 0.50000)
 adjoint: exp(ad(itx))Y
(-0.28033-0.50000)(-0.35355 0.60355)( 0.42678 0.00000)
( 0.35355-0.60355)( 0.14645 0.00000)(-0.35355-0.60355)
( 0.42678 0.00000)( 0.35355 0.60355)(-0.28033 0.50000)

 t =   3.141593
 Adjoint: Ad(itx)Y = exp(-itx) Y exp(itx)
( 0.00000-0.00000)( 0.00000-0.00000)( 1.00000 0.00000)
(-0.00000 0.00000)(-1.00000 0.00000)( 0.00000 0.00000)
( 1.00000 0.00000)(-0.00000-0.00000)( 0.00000 0.00000)
 adjoint: exp(ad(itx))Y
(-0.00000 0.00000)( 0.00000-0.00000)( 1.00000 0.00000)
(-0.00000 0.00000)(-1.00000 0.00000)( 0.00000 0.00000)
( 1.00000 0.00000)(-0.00000-0.00000)(-0.00000-0.00000)

続行するには何かキーを押してください . . .