バネ埋め込み法を使って隣接行列からグラフを描画する

5 min read読了の目安(約4500字

"グラフからコミュニティ構造を抽出する 〜リッチフローによるグラフの時間発展〜" ではリッチフローによる隣接行列を時間変化を計算するプログラムを実装しました。

しかし隣接行列のリストが得られたとしても、そこから直ちに上記のような可視化を行うことはできません。この記事では隣接行列からどうやってグラフの可視化を行うかについて解説していきたいと思います。

ベクトルや行列の実装は前の記事で実装したものを使います。この記事で必要になる実装は以下のトグルに記載しておきました。

線形代数の実装
import Data.Vector.Sized (Vector)
import qualified Data.Vector.Sized as V

--------------
-- ベクトル
--------------

-- | ベクトルのスカラー倍
(*^) :: Num a => a -> Vector n a -> Vector n a
a *^ v = fmap (a*) v

-- | ベクトル間の足し算
(^+^) :: Num a => Vector n a -> Vector n a -> Vector n a
(^+^) = V.zipWith (+)

-- | ベクトル間の引き算
(^-^) :: Num a => Vector n a -> Vector n a -> Vector n a
(^-^) = V.zipWith (-)

-- | ベクトルのL2ノルム
norm2 :: Floating a => Vector n a -> a
norm2 = sqrt . V.sum . fmap (^2)

-- | ベクトル間のユークリッド距離
distance :: Floating a => Vector n a -> Vector n a -> a
distance xs ys = norm2 (xs ^-^ ys)

--------------
-- 行列
--------------

-- | m × n 行列
type Matrix m n a = Vector m (Vector n a)

-- | 全て同じ値を持つ行列を作成する関数
konst :: (KnownNat m, KnownNat n) => a -> Matrix m n a
konst a = V.replicate (V.replicate a)

-- 行列のスカラー倍
(*!!) :: Semiring a => a -> Matrix m n a -> Matrix m n a
a *!! m = mmap (otimes a) m

-- | 各成分ごとの演算
elementWise :: (a -> b -> c) -> Matrix m n a -> Matrix m n b -> Matrix m n c
elementWise op = V.zipWith (V.zipWith op)

-- | 行列の和
(!+!) :: Semiring a => Matrix m n a -> Matrix m n a -> Matrix m n a
(!+!) = elementWise oplus

バネ埋め込み法

バネ埋め込み法はグラフのノードを重り、エッジをバネだと思って、物理シミュレーションを行って安定した配置を描画するというグラフ描画手法の一つです[1]

ここではグラフのノードの重さは全て1であるとし、ノードには以下のような力が働くと仮定します。

バネの力

ノードiとノードjの間に長さlのエッジが存在した時、ノードiには以下の大きさを持つバネの力がノードjの方向に向かって働きます。

-k(l - d_{ij})

ここでkはバネ定数であり、d_{ij}はノードi, j間の距離です。

この力によって各ノードはエッジの重さが距離に対応するような位置に配置されることになります。つまり重いエッジを共有するノードはどんどん離れていき軽いエッジを共有するノードはどんどん近づくことになります。

クーロン力

バネの力だけだとノードとノードが重なって見えづらくなってしまう可能性があります。なのでノード間には以下の大きさを持つクーロン力(斥力)が働くと仮定します。

\frac{q}{d_{ij}^2}

qはクーロン力の大きさを表す定数です。この斥力は距離の二乗に反比例して小さくなっていきます。

摩擦力

上記2つの力だけだと振動して永遠に収束しない可能性があるのでノードには摩擦力がかかり徐々に減衰していくとします。

\mu

\muは摩擦係数で常にノードが動く方向と反対方向に掛かる力です。なので単純にノードの速度を減速させていくような作用と考えることができます。

原点への引力

バネ埋め込み法で調べると上記3つの力が基本的なものとして説明されていますが、これだけだとグラフが連結でなかった場合に非連結成分同士が無限に離れていってしまうので原点に対して少し引力が働くようにしておくとグラフの配置が原点を中心とするようになり計算が安定します。

-\frac{g}{d_{i0}^2}

gは原点への引力の係数であり、d_{i0}は原点からノードまでの距離です。


以上4つの力がグラフのノードに働くとした上で適当な初期値から始めて、時間変化が予め決めた閾値より低くなるまで、オイラー法で時間発展を計算すればグラフの配置が求まるという寸法です。この考え方はグラフを埋め込む空間の次元には依存しないので平面でも3次元空間に埋め込む場合でも同様の方法を用いることができます。

実装すると以下のようなコードになります。

-- | ベクトルのノルムを1にする
normalize' :: Floating a => Vector n a -> Vector n a
normalize' v = fmap (/ norm2 v) v

-- | バネ埋め込み法によるグラフ描画
springEmbed :: (KnownNat m, KnownNat n)
            => Matrix m m Double  -- 隣接行列
            -> Matrix m n Double  -- 初期位置
            -> Matrix m n Double  -- 計算位置
springEmbed adj ps0 = flip fix (ps0, konst 0.0) \loop (ps, vs) ->
  let vs' = flip V.imap vs \i v ->
          let p = ps `V.index` i
              f j p' accum = if i == j then accum
                  else let w  = adj `V.index` i `V.index` j
                           d  = distance p p'
                           f1 = q / d^2                                -- クーロン力
                           f2 = if w == 0.0 then 0.0 else -k * (d - w) -- バネの力
                           f3 = -g / d^2                               -- 原点への引力
                        in accum ^+^ ((f1 + f2) *^ normalize' (p ^-^ p')) ^+^ (f3 *^ normalize' p)
              dv = V.ifoldr f (V.replicate 0) ps
           in mu *^ (v ^+^ (dt *^ dv))
      ps' = ps !+! (dt *!! vs')
   in if norm2 (fmap norm2 vs') < e then ps' else loop (ps', vs')
  where
    g  = 0.01 -- 原点への引力
    q  = 0.1  -- ノード間の斥力
    k  = 1.0  -- ばね定数
    mu = 0.9  -- 摩擦係数
    dt = 0.1  -- 微小時間
    e  = 0.1  -- 収束判定のしきい値

計算結果 Matrix m n am個のノードのn次元空間での座標がまとめられた行列です。最初は適当な初期位置から計算を始めれば良いですし、もし隣接行列の時間発展を扱っている場合は、一つ前の時間における配置を初期値にすることで滑らかなグラフの時間発展を描画することができるでしょう。あとはgloss等を用いて求めた配置にノードとエッジを描画すればグラフの可視化の完成です 👏

P.S. 参考までにglossの使い方を載せた記事へのリンクを張っておきます

https://qiita.com/lotz/items/eb73e62a64bc208c2dd6

\読んでいただきありがとうございました!/
この記事が面白かったら いいね♡ をいただけると嬉しいです☺️
100円からでも サポート¥ をいただければ次の記事を書くため励みになります🙌

脚注
  1. 理論的な話はこちら Tutte embedding - Wikipedia ↩︎

この記事に贈られたバッジ