行列とベクトルに関する計算: Sample 1: 行列と行列の掛け算(FORTRAN)

言語の変更： C版

■ 概要

行列と行列の掛け算は、特殊な場合を除いて、計算量のオーダーがO(N³)と決まっています。ここではある特殊な場合を考えます。つまり行列と行列を掛け算したあとにその結果の行列はまたベクトルを掛け算したい場合です。

式で書くと式(1)になります。この計算に必要なのはベクトルyだけであり、行列Cはyを計算するためだけであり、yが求まったらCは不要とします。

... 式(1)

式(1)を素直にコーディングしたのはCode 1になります。式(1)の前半はオーダーがO(N³)の行列と行列の掛け算、後半はオーダーがO(N²)の行列とベクトルの掛け算。全体のオーダーは最大のオーダーで決まるので、O(N³)。

さて式(1)をまとめて書くと、y=(AB)xとなりますが、行列・ベクトルの掛け算の順番を入れ替えることはできるので、式(1)はy=A(Bx)とも書けます。さらに2つの式に分けると式(2)になります。

... 式(2)

式(2)でご注意していただきたいのは行列と行列の掛け算がなくなることです。その代わりに2つの行列とベクトルの掛け算がありますが、行列とベクトルの掛け算はオーダーがO(N²)なので、式(2)のオーダーは式(1)より一つ低く、O(N²)です。

この次にCode 1とCode 2を実際に実行するときにどのぐらい計算時間の差が出るかを見てみましょう。

■ ソースコード

◆ Code 1 ◆ Code 2

 
      program main
      implicit none
      integer i, j, k, n, m
      real*8, allocatable :: a(:,:), b(:,:), c(:,:), x(:), y(:)
      integer*8 time0, time1, dtime
      real*8 time
c
      write(*,*)"Matrix size    Elapsed time [sec]"
      do m=1,20
c Matrix size
        n = m * 100

c Allocation
        allocate(a(n,n), b(n,n), c(n,n), x(n), y(n))

c Initialization
        do i=1,n
          do j=1,n
            a(j,i) = 1d0
            b(j,i) = 2d0
          end do
          x(i) = 1d0
        end do

c Start time
        call system_clock(time0)

c Main calculation 1: matrix-matrix product
        do i=1,n
          do j=1,n
            c(j,i) = 0d0
            do k=1,n
              c(j,i) = c(j,i) + a(j,k) * b(k,i)
            end do
          end do
        end do

c Main calculation 2: matrix-vector product
        do i=1,n
          y(i) = 0d0
          do k=1,n
            y(i) = y(i) + c(k,i) * x(k)
          end do
        end do

c Finish time
        call system_clock(time1, dtime)

c Output time
        time = 1d0*(time1-time0)/dtime
        write(*,"(i12,f16.7)")n, time

c Deallocation
        deallocate(a, b, c, x, y)
      end do

      end program

 
      program main
      implicit none
      integer i, j, k, n, m
      real*8, allocatable :: a(:,:), b(:,:), x(:), y(:), z(:)
      integer*8 time0, time1, dtime
      real*8 time
c
      write(*,*)"Matrix size    Elapsed time [sec]"
      do m=1,15
c Matrix size
        n = m * 1000

c Allocation
        allocate(a(n,n), b(n,n), x(n), y(n), z(n))

c Initialization
        do i=1,n
          do j=1,n
            a(j,i) = 1d0
            b(j,i) = 2d0
          end do
          x(i) = 1d0
        end do

c Start time
        call system_clock(time0)

c Main calculation 1: matrix-vector product
        do i=1,n
          z(i) = 0d0
          do k=1,n
            z(i) = z(i) + b(k,i) * x(k)
          end do
        end do



c Main calculation 2: matrix-vector product
        do i=1,n
          y(i) = 0d0
          do k=1,n
            y(i) = y(i) + a(k,i) * z(k)
          end do
        end do

c Finish time
        call system_clock(time1, dtime)

c Output time
        time = 1d0*(time1-time0)/dtime
        write(*,"(i12,f16.7)")n, time

c Deallocation
        deallocate(a, b, x, y, z)
      end do

      end program

■ 計算時間の測定結果

Code 1とCode 2を実行したときの計算時間をそれぞれ青線と赤線で図1に示します。

Code 2の計算時間はCode1に比べて非常に短いので、Code 1と同じサイズNで測定すると、計算時間が測定できません。ここではCode 2のサイズNをCode 1より大きい範囲にしています。

図1 測定時間

■ 考察

図1より解かりますが、Code 1の計算時間はCode 2より非常に長いです。しかもサイズNが大きくなるにつれて、その差がどんどん開きます。ちなみにN=2000ではCode 1の計算時間は既にCode 2の5000倍程度になっています。Code 1をNが10000まで計算してもよいですが、1日ぐらいかかりそうなので、止めました。この結果より、O(N³)の恐ろしさを実感できるのではないでしょうか。

前: 行列とベクトルに関する計算

次: ダウンロード

行列とベクトルに関する計算: Sample 1: 行列と行列の掛け算(FORTRAN)