Fruchterman, Reingold のグラフ描画アルゴリズム

はしがき

1991 年に発表されたグラフ描画アルゴリズムです。
こちらが掲載されている論文も今日現在で Google Scholar 調べの引用数が 1592 と凄いです。

アルゴリズム

このアルゴリズムも Kamada, Kawai の方法と同様に、最初に頂点の位置を決定し、その後、頂点同士を結ぶ直線で辺を表す形のアルゴリズムです。

頂点の位置は、次の二つの原理に従うように決定されます。

  • 頂点同士が近過ぎない。
  • 隣接する頂点がそれなりに近い。

より具体的には、各頂点の初期位置をランダムに決定した後、アルゴリズムはある回数だけ次の過程を繰り返します。

  1. 各頂点に働く斥力をヒューリスティックに計算する。
  2. 各辺の端点同士に働く引力をヒューリスティックに計算する。
  3. 引力、斥力、温度を元に各頂点の次の位置をヒューリスティックに計算する。
  4. 描画領域からはみ出した頂点の位置を修正する。
  5. 温度を下げる。

各ステップは入力グラフの頂点数の二乗のオーダの時間計算量で計算することができます。

後でプログラムのサンプルを書く

Kamada, Kawai のグラフ描画アルゴリズム

はしがき

今回はフロイド法のところで言及したグラフ描画アルゴリズムについて書きます。
私はグラフ描画の世界に関しては何も知らない素人で、このアルゴリズムも古典的な有名どころらしきものを適当に選んだだけなのですが、このアルゴリズムが掲載されている An Algorithm for Drawing General Undirected Graphs は今日の時点で Google Scholar 調べで 1198 の記事から引用されているというとても凄い論文です。

アルゴリズム

アルゴリズムはとてもシンプルなもので、グラフの各頂点の現在の描画位置を引数にとるエネルギー関数を定義し、その極小点を Newton-Raphson 法で求めるだけです。
頂点を結ぶ辺は、頂点の位置が定まった後に直線によって描画されます。
よって（入力が平面グラフであっても）辺は交差している場合があります。

エネルギー関数の詳細と擬似コードは後で書く

C＋＋プログラムの一例です。
非常に長くなってしまいましたが、下手にソースを分割すると後で使いにくくなるおそれがあるので、過去に似たようなコードが出てきたような部分も含めて、全て一つのファイルに収めました。
また計算結果を出力をどうするか迷ったのですが、HTMLとJavaScriptのコードを出力させてブラウザに表示するようにしました。
一般的なブラウザさえあれば、余計なソフトを追加でインストールしなくても大丈夫です。

ちなみに、このアルゴリズムは定義より明らかに、非連結なグラフに対して適用すると挙動がおかしくなります。
しかしランダムグラフで作っている入力グラフは平気で非連結になるので、サンプルの出力がおかしい場合は何回かリロードしてみてください。

#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iomanip>
#include <iostream>
using namespace std;

inline double roll() { return (double)rand() / RAND_MAX; }
const double MAX = pow(2.0, 64.0);
double** allocate_square_matrix(int num_element);

int main() {
	srand((unsigned int)time(0));
	
	const int num_vertices = 15;
	const int& N = num_vertices;
	
	double** adjacent = allocate_square_matrix(N);
	double**& A = adjacent;
	
	// make up an input graph (make a random graph)
	const double p = 0.2;
	const double max_length = 20;
	const double min_length = 1;
	for (int v = 0; v < N; v++) {
		A[v][v] = MAX;
		for (int u = v+1; u < N; u++) {
			double l = (roll() < p)
			? (max_length - min_length) * roll() + min_length : MAX;
			A[v][u] = l;
			A[u][v] = l;
		}
	}
	
	cerr << "The adjacent matrix of the input graph" << endl;
	cerr.setf(ios::fixed, ios::floatfield);
	cerr << setprecision(1);
	for (int v = 0; v < N; v++) {
		for (int u = 0; u < N; u++) {
			cerr << setw(4);
			if (A[v][u] == MAX)
				cerr << " MAX ";
			else
				cerr << A[v][u] << " ";
		}
		cerr << endl;
	}
	cerr << endl;

  // floyd's method
  double** distance = allocate_square_matrix(N);
  double**& D = distance;
  for (int v = 0; v < N; v++)
    for (int u = 0; u < N; u++)
      if (v == u)
        D[v][u] = 0.0;
      else
        D[v][u] = A[v][u];

  for (int k = 0; k < N; k++)
    for (int v = 0; v < N; v++)
      for (int u = 0; u < N; u++)
        if (D[v][k] + D[k][u] < D[v][u])
          D[v][u] = D[v][k] + D[k][u];

	cerr << "The distance matrix" << endl;
	cerr.setf(ios::fixed, ios::floatfield);
	cerr << setprecision(1);
	for (int v = 0; v < N; v++) {
		for (int u = 0; u < N; u++) {
			cerr << setw(4);
			if (D[v][u] == MAX)
				cerr << " MAX ";
			else
				cerr << D[v][u] << " ";
		}
		cerr << endl;
	}
	cerr << endl;

  // Kamada and Kawai's algorithm
  double** l = allocate_square_matrix(N);
  double** k = allocate_square_matrix(N);
  double L0 = 500;
  double K = 10;
  double max_distance = 0;
  for (int v = 0; v < N; v++)
    for (int u = 0; u < N; u++)
      if (D[v][u] > max_distance)
        max_distance = D[v][u];
  double L = L0 / max_distance;
  for (int v = 0; v < N; v++)
    for (int u = 0; u < N; u++) {
      l[v][u] = L * D[v][u];
      k[v][u] = K / D[v][u] / D[v][u];
    }

  double* x = new double[N];
  double* y = new double[N];
  for (int v = 0; v < N; v++) {
    x[v] = L0 * roll();
    y[v] = L0 * roll();
  }

	cerr.setf(ios::scientific, ios::floatfield);
	cerr << setprecision(1);
  const double epsilon = 0.001;
  const double& E = epsilon;
  while (1) {
    double max_delta = 0;
    double max_dEx = 0;
    double max_dEy = 0;
    int argmax_delta = -1;
    for (int v = 0; v < N; v++) {
      double dEx = 0;
      double dEy = 0;
      for (int u = 0; u < N; u++)
        if (u != v) {
          double d = (x[v] - x[u]) * (x[v] - x[u]) + 
            (y[v] - y[u]) * (y[v] - y[u]);
          d = sqrt(d);
          dEx += k[v][u] * ((x[v] - x[u]) - l[v][u] * (x[v] - x[u]) / d);
          dEy += k[v][u] * ((y[v] - y[u]) - l[v][u] * (y[v] - y[u]) / d);
        }
      double delta = sqrt(dEx * dEx + dEy * dEy);
      if (delta > max_delta) {
        max_delta = delta;
        max_dEx = dEx;
        max_dEy = dEy;
        argmax_delta = v;
      }
    }
    if (max_delta <= epsilon)
      break;
    //cerr << max_delta << endl;

    double delta = max_delta;
    double dEx = max_dEx;
    double dEy = max_dEy;
    int v = argmax_delta;
    while (delta > epsilon) {
      double dEx2 = 0;
      double dEy2 = 0;
      double dExy = 0;
      double dEyx = 0;
      for (int u = 0; u < N; u++)
        if (u != v) {
          double d = (x[v] - x[u]) * (x[v] - x[u]) + 
            (y[v] - y[u]) * (y[v] - y[u]);
          d *= sqrt(d);
          dEx2 += k[v][u] * (1 - l[v][u] * (y[v] - y[u]) * (y[v] - y[u]) / d);
          dExy += k[v][u] * l[v][u] * (y[v] - y[u]) * (x[v] - x[u]) / d;
          dEyx += k[v][u] * l[v][u] * (x[v] - x[u]) * (y[v] - y[u]) / d;
          dEy2 += k[v][u] * (1 - l[v][u] * (x[v] - x[u]) * (x[v] - x[u]) / d);
        }

      double a, b, c, d, det;
      det = dEx2 * dEy2 - dExy * dEyx;
      a = dEy2 / det;
      b = -dExy / det;
      c = -dEyx / det;
      d = dEx2 / det;
      x[v] += -dEx * a + -dEy * b;
      y[v] += -dEx * c + -dEy * d;

      dEx = 0;
      dEy = 0;
      for (int u = 0; u < N; u++)
        if (u != v) {
          double d = (x[v] - x[u]) * (x[v] - x[u]) + 
            (y[v] - y[u]) * (y[v] - y[u]);
          d = sqrt(d);
          dEx += k[v][u] * ((x[v] - x[u]) - l[v][u] * (x[v] - x[u]) / d);
          dEy += k[v][u] * ((y[v] - y[u]) - l[v][u] * (y[v] - y[u]) / d);
        }
      delta = sqrt(dEx * dEx + dEy * dEy);
      //cerr << "  " << delta << endl;
    }
  }

  // HACK
  for (int v = 0; v < N; v++) {
    x[v] += 50;
    y[v] += 50;
  }

	cout.setf(ios::fixed, ios::floatfield);
	cout << setprecision(3);
  cout << "<html>" << endl;
  cout << "  <head>" << endl;
  cout << "    <title>A graph</title>" << endl;
  cout << "    <script type='text/javascript'>" << endl;
  cout << "    onload = function() {" << endl;
  cout << "      var canvas = document.getElementById('graphcanvas');" << endl;
  cout << "      if (!canvas || !canvas.getContext)" << endl;
  cout << "        return false;" << endl;
  cout << "      var ctx = canvas.getContext('2d');" << endl;
  cout << "      ctx.fillStyle = '#ffffff';" << endl;
  cout << "      ctx.strokeStyle = '#000000';" << endl;
  cout << "      ctx.textAlign = 'center';" << endl;
  for (int v = 0; v < N; v++)
    for (int u = v+1; u < N; u++)
      if (A[v][u] < MAX) {
        cout << "      ctx.beginPath();" << endl;
        cout << "      ctx.moveTo(" << x[v] << ", " << y[v] << ");" << endl;
        cout << "      ctx.lineTo(" << x[u] << ", " << y[u] << ");" << endl;
        cout << "      ctx.closePath();" << endl;
        cout << "      ctx.stroke();" << endl;
      }
  for (int v = 0; v < N; v++) {
    cout << "      ctx.beginPath();" << endl;
    cout << "      ctx.arc(" << x[v] << ", " << y[v] 
         << ", 10, 0, Math.PI*2, true);" << endl;
    cout << "      ctx.fill();" << endl;
    cout << "      ctx.stroke();" << endl;
    cout << "      ctx.strokeText('" << v
         << "', " << x[v] << ", " << (y[v]+4) << ");" << endl; // HACK
  }
  cout << "    };" << endl;
  cout << "    </script>" << endl;
  cout << "  </head>" << endl;
  cout << "  <body>" << endl;
  cout << "    <canvas id='graphcanvas' width='600' height='600'></canvas>" << endl;
  cout << "  </body>" << endl;
  cout << "</html>" << endl;

  return 0;
}

double** allocate_square_matrix(int num_element) {
	double** matrix = new double*[num_element];
	matrix[0] = new double[num_element * num_element];
	for (int v = 1; v < num_element; v++)
		matrix[v] = matrix[v-1] + num_element;
  return matrix;
}

フロイド法

はしがき

先日、某有名グラフ可視化ソフトがあまりにもオモチロイ図を吐き出してくださるので、自分で適当な可視化アルゴリズムを勉強して実装してみようと思い立ちました。
そこで適当な論文を選んで読んでいたところ、サブルーチンとしてフロイド法が出てきたので本日はこれについてまとめてみようと思いました。
そのうちそっちのグラフ可視化アルゴリズムについても書くつもりです。

概要と歴史

フロイド法はグラフの全点対最短経路問題を解くためのアルゴリズムです。
グラフは有向／無向のどちらでもよいです。
グラフが負の長さのサイクルを持たないとき、このアルゴリズムは全ての頂点間の最短経路の長さを保持する配列を作成します*1。
グラフが負の長さのサイクルを持つ場合にはそのサイクルを通過する経路を持つ頂点間の最短経路が正しく定義できなくなってしまうのですが、上記の配列のある要素を調べると、グラフが負の長さのサイクルを持つ場合にはそれを検出することができます。
グラフが負の長さのサイクルを持つ場合、負の長さのサイクルとは無関係な頂点間の最短経路がどうなるかは私はまだきちんと考えてはいません。

フロイド法はワーシャル・フロイド法やロイ・ワーシャル法などとも呼ばれるらしいです。
ワーシャル（Stephan Warshall）は 1962 年 1 月に出版された A Theorem on Boolean Matrices でグラフの連結性を計算するアルゴリズムを提案しました。
このアルゴリズムは表題通り二値の（隣接）行列を入力に取るのですが、アルゴリズムの構造を保ったまま、入力を辺の長さの隣接行列に置き換え、中の演算を距離の計算に置き換えたものがフロイド法です。
1962 年 6 月に出版されたフロイド（Robert Floyd）の Algorithm97 Shortest Path では上記の論文がきちんと引用されています。
1959 年に出版されたロイ（Bernard Roy）の論文はフランス語なので私には読めませんでしたが、タイトルから何となくワーシャルと同様に連結性に関するアルゴリズムなのではないかと思います。

アイディアと擬似コード

グラフの頂点数を n とし、各頂点に 1 から n までの番号を割り振ります。
全ての頂点対 i, j に対して、i, j 以外には頂点 1, 2, ..., k 以外の頂点を通らない経路のうち長さが最小のものを計算し、それを利用して次は、全ての頂点対 i, j に対して、i, j 以外には頂点 1, 2, ..., k, k+1 以外を通らない経路のうち長さが最小のものを計算する、という処理を k = 0 から n まで順に行うのがフロイド法です。
k = 0 のときは、入力のグラフにおいて i, j 間に辺が張られているかどうかを調べます。
k = K > 0 のときは、i, j 間の最短経路は、k = K-1 のときの i, j 間の最短経路と、k = K-1 における i, k 間の最短経路と k, j 間の最短経路を合わせたもののうち、長さが短いほうです。
フロイド法の時間計算量は O(n^3) となります。

後で擬似コードを書く

C＋＋プログラムの一例です。
入力グラフはダイクストラ法のプログラム例と同様にランダムグラフを用いて適当に作りました。

#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iomanip>
#include <iostream>
using namespace std;

inline double roll() { return (double)rand() / RAND_MAX; }
const double MAX = pow(2.0, 64.0);
double** allocate_square_matrix(int num_element);

int main() {
	srand((unsigned int)time(0));
	
	const int num_vertices = 5;
	const int& N = num_vertices;
	
	double** adjacent = allocate_square_matrix(N);
	double**& A = adjacent;
	
	// make up an input graph (make a random graph)
	const double p = 0.3;
	const double max_length = 20;
	const double min_length = 1;
	for (int v = 0; v < N; v++) {
		A[v][v] = MAX;
		for (int u = v+1; u < N; u++) {
			double l = (roll() < p)
			? (max_length - min_length) * roll() 
			+ min_length : MAX;
			A[v][u] = l;
			A[u][v] = l;
		}
	}
	
	cout << "The adjacent matrix of the input graph" << endl;
	cout.setf(ios::fixed, ios::floatfield);
	cout << setprecision(1);
	for (int v = 0; v < N; v++) {
		for (int u = 0; u < N; u++) {
			cout << setw(4);
			if (A[v][u] == MAX)
				cout << " MAX ";
			else
				cout << A[v][u] << " ";
		}
		cout << endl;
	}
	cout << endl;

  // floyd's method
  double** distance = allocate_square_matrix(N);
  double**& D = distance;
  for (int v = 0; v < N; v++)
    for (int u = 0; u < N; u++)
      if (v == u)
        D[v][u] = 0.0;
      else
        D[v][u] = A[v][u];

  for (int k = 0; k < N; k++)
    for (int v = 0; v < N; v++)
      for (int u = 0; u < N; u++)
        if (D[v][k] + D[k][u] < D[v][u])
          D[v][u] = D[v][k] + D[k][u];

	cout << "The distance matrix" << endl;
	cout.setf(ios::fixed, ios::floatfield);
	cout << setprecision(1);
	for (int v = 0; v < N; v++) {
		for (int u = 0; u < N; u++) {
			cout << setw(4);
			if (D[v][u] == MAX)
				cout << " MAX ";
			else
				cout << D[v][u] << " ";
		}
		cout << endl;
	}
	cout << endl;

  return 0;
}

double** allocate_square_matrix(int num_element) {
	double** matrix = new double*[num_element];
	matrix[0] = new double[num_element * num_element];
	for (int v = 1; v < num_element; v++)
		matrix[v] = matrix[v-1] + num_element;
  return matrix;
}

*1:この手のアルゴリズムで、アルゴリズムの停止後、O(1) の時間計算量で経路が得られるように書いている記事は、厳密には間違いである場合が多いと思います。O(1) で経路が得られるようにアルゴリズムを改良するのは大抵の場合は簡単なのですが、その多くの場合は時間計算量と空間計算量が元のアルゴリズムより増加してしまいます。フロイド法の場合は、倍のメモリを使うことで、アルゴリズムの停止後 O(n) である頂点間の最短経路を出力するよう、時間計算量を変えずにアルゴリズムを改良することはできます。

ダイクストラ法

はしがき

このブログは、ウィキペディアの日本語ページは無い程度にはマイナーな知識を書き記しておくために始めたものですが、最初のうちは練習がてら、既に広く知られていることについても書いてみようと思います。

概要と歴史

ダイクストラ法は負の長さを持たないグラフ（有向／無向問わず）の最短経路問題を解くためのアルゴリズムです。

元々 1959 年にエドガー・ダイクストラ（Edsger Dijkstra）が発表したもの*1はグラフ上の二頂点間の最短経路を求めるためのアルゴリズムだったのですが、現在では、二頂点間の最短経路の長さを求めるアルゴリズムや、ある頂点からの最短経路木を構成するアルゴリズムをダイクストラ法と呼んでいる資料も数多く存在します（なぜならこれらは全て、ダイクストラが提案したものとほぼ同じ方法で求めることができるからです）。
また、ダイクストラが発表したアルゴリズムは、近年の研究であれば当たり前に求められるレベルの詳細すらも示されていない単なる「アイディア」なのですが、その詳細を何らかの方法で定めた具体的なアルゴリズムもまとめてダイクストラ法と呼ばれることがあるようです。

アイディアと擬似コード・プログラムの例

始点からある頂点 v への最短経路は、始点から v と隣接しているある頂点 u への最短経路に辺 { u, v } を付け足したものである、という性質がグラフ上の始点以外の全ての頂点について成立します（始点から始点への最短経路は空経路とします）。
この性質を利用して、始点からそれぞれの頂点への最短経路の長さを、始点に近いほうから順に、連鎖的に定めていこう、というのがダイクストラ法のアイディアです。
終点への最短経路の長さが求められたら、上の性質を利用して、今度は逆向きに頂点をたどっていけば最短経路自体を得ることができます。

あとで擬似コードを書く

擬似コードよりも先にC＋＋プログラムの一例を書いてしまいました。
入力グラフは適当にランダムグラフを作って使っています。

#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iomanip>
#include <iostream>
using namespace std;

inline double roll() { return (double)rand() / RAND_MAX; }
const double MAX = pow(2.0, 64.0);

int main() {
  srand((unsigned int)time(0));

  const int num_vertices = 10;
  const int& N = num_vertices;

  double** adjacent = new double*[N];
  double**& A = adjacent;
  A[0] = new double[N * N];
  for (int v = 1; v < N; v++)
    A[v] = A[v-1] + N;

  // make up an input graph (make a random graph)
  const double p = 0.3;
  const double max_length = 20;
  const double min_length = 1;
  for (int v = 0; v < N; v++) {
    A[v][v] = MAX;
    for (int u = v+1; u < N; u++) {
      double l = (roll() < p)
        ? (max_length - min_length) * roll() + min_length : MAX;
      A[v][u] = l;
      A[u][v] = l;
    }
  }

  cout << "The adjacent matrix of the input graph" << endl;
  cout.setf(ios::fixed, ios::floatfield);
  cout << setprecision(1);
  for (int v = 0; v < N; v++) {
    for (int u = 0; u < N; u++) {
      cout << setw(4);
      if (A[v][u] == MAX)
        cout << " MAX ";
      else
        cout << A[v][u] << " ";
    }
    cout << endl;
  }
  cout << endl;

  // dijkstra's method

  // initialize
  double* length = new double[N];
  double*& L = length;
  bool* done = new bool[N];
  bool*& D = done;
  int* parent = new int[N];
  int*& P = parent;
  for (int v = 0; v < N; v++) {
    L[v] = MAX;
    D[v] = false;
    P[v] = -1;
  }

  // make shortest path tree
  L[0] = 0;
  P[0] = 0;
  while (1) {
    int min = -1;
    for (int v = 0; v < N; v++)
      if (!D[v] && (min == -1 || L[v] < L[min]))
        min = v;
    if (min == -1)
      break;

    D[min] = true;
    for (int v = 0; v < N; v++)
      if (L[min] + A[min][v] < L[v]) {
        L[v] = L[min] + A[min][v];
        P[v] = min;
      }
  }

  cout.setf(ios::fixed, ios::floatfield);
  cout << setprecision(1);
  for (int v = 0; v < N; v++) {
    cout << "L[" << v << "] = " << setw(5);
    if (L[v] == MAX)
      cout << "MAX" << endl;
    else
      cout << L[v] << ", P[" << v << "] = " << P[v] << endl;
  }
  cout << endl;

  // output a path
  int curr = rand() % N;
  if (D[curr]) {
    cout << "A path from " << curr << " to 0" << endl;
    while (curr != 0) {
      cout << curr << ", ";
      curr = P[curr];
    }
    cout << 0 << endl;
  }

  return 0;
}

書き方がC＋＋らしくないのはただの趣味です。

実装のちがいによる時間計算量の変化

ダイクストラ法は

  • 入力グラフをどう与えるか
  • 始点からの経路の長さがまだ確定していない頂点をどう探すか

によって時間計算量が変化します。

まず一般に、ダイクストラ法の最短経路木を作る部分の時間計算量は O(I + \sum (S_i + U_i) です。
ただし I は経路長用配列の初期化、S_i と U_i はそれぞれ i 回目のループにおいて未確定頂点の探索と経路長の更新にかかる時間計算量です。
グラフの頂点数を n、辺数を m とすると、 I = O(n) であり、ループの回数は最大で n なので（最短経路木を構成する場合はもちろん、終点までの最短経路が発見された直後に打ち切る場合でも最短経路が全ての頂点を通過する可能性があります） \sum (S_i + U_i) は I よりも大きくなります。

入力グラフが隣接グラフによって与えられるとき、 U_i は常に \Omega(n) です。
また、 S_i は高々 O(n) （力任せに探す場合）です。
よって入力グラフが隣接グラフによって与えられる場合の時間計算量は O(I + \sum (S_i + U_i)) = O(n^2) となります。
上記のプログラムはこの形で実装してあります。

入力が隣接リストによって与えられる場合には、 O(\sum U_i) = O(m L) です。
ただし L は未探索頂点を効率よく探索するためのデータ構造を管理するのにかかる時間計算量の上界です。
力任せに探す場合は L = O(1) なので、上と同様に O(I + \sum (S_i + U_i)) = O(n^2) です（ちなみに m \leq n^2 です。完全グラフの場合）。
頂点を最小ヒープ（二分ヒープ）によって管理する場合は、 S = O(\log n), L = O(\log n) なので、 O(I + \sum (S_i + U_i)) = O(n \log n + m \log n) となります。
頂点をフィボナッチヒープ*2や Brodal Queue*3などを用いると時間計算量を更に落とすことができます。

また、最短経路木を作ったあと、ある頂点から始点への経路を求める部分の時間計算量は明らかに O(n) となります。
ただし、最短経路木において各頂点の親がどの頂点かを保持しない場合には、入力が隣接行列の場合は O(n^2)、隣接リストの場合は O(m) となります。

上記のいかなる方法でも、入力グラフを表す部分を除いて、ダイクストラ法自体にかかる空間計算量は全て O(n) です。

今度 tex の部分を別のサービスに置き換えるかも

*1:A Note on Two Problems in Connexion with Graphs。ちなみに An Interview with Edsger W. Dijkstra によれば作ったのは発表の三年前の 1956 年だったそうです

*2:昔勉強したけど忘れました

*3:これから勉強する予定です