疎行列の計算を実装してグラフ理論をかじってみる

経路数を計算する

　それではおためし。単位行列Aとfig-2に示した隣接行列Bを用意し、A = A * Bを三回繰り返してB^3を求めます。

list-4　01_adjacency\01_adjacency.cpp

#include "sparse/matrix.h"
#include "sparse/multiplies.h"
#include "sparse/fill_diagonal.h"

#include <iostream>   // cout, endl

int main() {
  using namespace sparse;

  matrix<float> A;
  fill_diagonal(A, 4, 1.0f); // 4x4の単位行列
  
  matrix<float> B; // 隣接行列
  B.insert(0, 0, 1.0f);
  B.insert(0, 1, 1.0f);
  B.insert(0, 2, 1.0f);
  B.insert(1, 2, 1.0f);
  B.insert(2, 3, 1.0f);
  B.insert(3, 1, 1.0f);
  B.finalize();

  // A = A * B を3回行う
  A = multiplies(A, B, 4, 4, 4);
  A = multiplies(A, B, 4, 4, 4);
  A = multiplies(A, B, 4, 4, 4);

  A.enumerate([](int row, int col, float val) {
                std::cout << '[' << row << ',' << col << "] " 
                          << val << std::endl;}); 
  std::cout << std::endl;
}

　隣接行列を3回かけたときの成分（i,j）が何を意味するかというと、「頂点iから頂点jへ3つの辺をたどって到達する経路の数」なんです。例えば「[0,3] 2」とありますから、島0から3本の橋を渡って島3にたどり着くルートは2通りです。確かに0→0→2→3と0→1→2→3の2通りですね。

経路の最短／最長距離を求める

　さて、ここから面白くなりますよ。

　K行L列の行列AとL行M列の行列Bの掛け算関数multipliesは、その成分（row,col）を以下のアルゴリズムで求めています。

val = 0;
for ( i = 0; i < L; ++i ) {
  val = val + A(row, i) * B(i, col);

　for-loopの中で乗算と加算が1つずつ現れてますね。この乗算と加算を引数で与えるmultipliesを用意します。

list-5　02_adjacency\02_adjacency.cpp

template<typename T, typename U, typename MulOp, typename AddOp>
matrix<T> multiplies(const matrix<T>& A, const matrix<U>& B, 
                     int k, int l, int m,
                     MulOp nmul, AddOp nadd) {
  matrix<T> out;
  // out[row, col] = Σ(A[row,i] * B[i,col]
  for ( int row = 0; row < k; ++row ) {
    for ( int col = 0; col < m; ++col ) {
      T s; T* sp = nullptr;
      T v; T* vp;
      T a; T* ap;
      U b; U* bp;
      for ( int i = 0; i < l; ++i ) {
        ap = A.at(row, i, &a); // a = A[row,i]
        bp = B.at(i, col, &b); // b = B[i,col]
        vp = nmul(ap, bp, &v); // v = a * b
        sp = nadd(sp, vp, &s); // s = s + v;
      }
      if ( sp != nullptr ) {
        out.insert(row, col, s);
      }
    }
  }
  out.finalize();
  return out;
}

　ここで与える乗算／加算関数MulOo／AddOpは、以下のプロトタイプに従うこととします。

T* operation(const T* a, const T* b, T* c)

　*aと*bを演算子の両辺とし、演算結果を*cにセットしてcを返します。ただし左辺／右辺が無効値であったときはa,bはnullptrですし、演算結果が無効値ならばnullptrを返します。

　ここで乗算と加算を以下のように定義します:

• 乗算：両辺の和を返す。ただしどちらかが無効値であるなら無効値を返す。
• 加算：両辺の「大きくない方」を返す。ただし一方が無効値なら他方を、両辺とも無効値なら無効値を返す。

　隣接行列にも手を加え、各成分を（1ではなく）橋の長さに置き換えます。

list-6　03_shortestpath\03_shortestpath.cpp

#include "sparse/matrix.h"
#include "sparse/multiplies.h"
#include "sparse/fill_diagonal.h"

// 乗算 <*>
//  0 <*> 0 = 0
//  x <*> 0 = 0
//  0 <*> x = 0
//  x <*> y = x+y
template<typename T> 
struct nullable_multiplies {
  T* operator()(const T* ap, const T* bp, T* result) const {
    if ( ap == nullptr ) return nullptr;
    if ( bp == nullptr ) return nullptr;
    *result = (*ap) + (*bp);
    return result;
  }
};

// 加算 <+>
//   0 <+> 0 = 0
//   0 <+> x = x
//   x <+> 0 = x
//   x <+> y = min(x,y)
template<typename T> 
struct nullable_min {
  T* operator()(const T* ap, const T* bp, T* result) const {
    if ( ap == nullptr && bp == nullptr ) return nullptr;
    if ( ap == nullptr ) *result = *bp;
    else *result = ( bp == nullptr ) ? *ap : ((*ap) < (*bp) ? (*ap) : (*bp));
    return result;
  }
};

// 加算 <+>
//   0 <+> 0 = 0
//   0 <+> x = x
//   x <+> 0 = x
//   x <+> y = max(x,y)
template<typename T> 
struct nullable_max {
  T* operator()(const T* ap, const T* bp, T* result) const {
    if ( ap == nullptr && bp == nullptr ) return nullptr;
    if ( ap == nullptr ) *result = *bp;
    else *result = ( bp == nullptr ) ? *ap : ((*ap) > (*bp) ? (*ap) : (*bp));
    return result;
  }
};

#include <iostream>   // cout, endl

int main() {
  using namespace sparse;

  matrix<float> A;
  fill_diagonal(A, 4, 0.0f);

  matrix<float> B;
  B.insert(0, 0, 0.8f);
  B.insert(0, 1, 1.5f);
  B.insert(0, 2, 1.8f);
  B.insert(1, 2, 1.1f);
  B.insert(2, 3, 2.1f);
  B.insert(3, 1, 2.2f);
  B.finalize();

  nullable_multiplies<float> mulop;
  nullable_min<float>        addop;

  A = multiplies(A, B, 4, 4, 4, mulop, addop);
  A = multiplies(A, B, 4, 4, 4, mulop, addop);
  A = multiplies(A, B, 4, 4, 4, mulop, addop);

  A.enumerate([](int row, int col, float val) {
                std::cout << '[' << row << ',' << col << "] " 
                          << val << std::endl;}); 
  std::cout << std::endl;
}

　この演算で得られる成分（i,j）は、「3本の橋をたどって島iから島jに到達するときの最短路長」なんです。加算の定義を「小さくない方」に差し替えれば最長路長が求まります。