VBAによる部分グラフ同型判定プログラムについて

Question

前提・実現したいこと

あるグラフの隣接行列A1とあるグラフの隣接行列A2がわかっているとして

A1 = P * A2__T__ * P__T__ (__T__は転置)

を満たす置換行列(各行各列に1が一つ、それ以外が0の行列)Pが存在するかどうか調べるプログラムを作成したいです。

置換行列Pが存在すれば、A2はA1に含まれる部分グラフだということまでわかったのですが、

Pの求め方がわからなくて困っています。

下のソースは置換行列Pが存在する場合の行列をA1とA2に設定しています。

ご教授願います。

よろしくお願いいたします。

該当のソースコード

Option Explicit

Sub main()
    Dim A1
    Dim A2

    A1 = [{0,1,1,1,1,0,0,0,0,0,0;1,0,0,0,0,1,0,0,0,0,0;1,0,0,0,0,0,1,1,1,0,0;1,0,0,0,0,0,0,0,0,0,0;1,0,0,0,0,0,0,0,0,0,0;0,1,0,0,0,0,0,0,0,2,1;0,0,1,0,0,0,0,0,0,0,0;0,0,1,0,0,0,0,0,0,0,0;0,0,1,0,0,0,0,0,0,0,0;0,0,0,0,0,2,0,0,0,0,0;0,0,0,0,0,1,0,0,0,0,0}]
    A2 = [{0,1,0,0,0;1,0,1,0,0;0,1,0,2,1;0,0,2,0,0;0,0,1,0,0}]

    ' 置換行列Pが存在すれば、Pありと出力
    If SubgraphHomotyping(A1, A2) = True Then
        Debug.Print "Pあり"
    Else
        Debug.Print "Pなし"
    End If
End Sub

' 置換行列Pが存在するか判定する処理(存在すればTrue, 存在しなければFalse)
Function SubgraphHomotyping(A1 As Variant, A2 As Variant) As Boolean

    ' 処理をここに記述

    If Pが存在する Then
        SubgraphHomotyping = True
    Else
        SubgraphHomotyping = False
    End If

End Function

Answer 1

http://www-erato.ist.hokudai.ac.jp/lecture2013/docs/washio2.pdf

こんなPDFを眺めてみました。p.7にある部分グラフ同型判定のところを見ますとNP完全問題とありました。p.9以降に様々な具体的手法が紹介されていますね。

p.9にUllmannの方法として「一定の約束事（行や列の順列など）で重複なく置換行列Pを生成する」と有りますが、この方法ではNP完全問題をまともに解く(有り得る様々なPを想定してしらみつぶしに調べてゆく)ということのように思えました。しかしこのページ以降に置換行列Pの大きさをしぼれそうな様々な手法が書かれていますね。

いずれにせよ少なくとも重複なく置換行列Pを生成するという論理が必要と思います。A1の大きさと重みありの個数を見ますと、まともに有り得る置換行列の全ての組み合わせを探索しても9!×2 = 725760となりますので素朴に総当たりしても(ご質問にある規模の問題なら)時間的には大丈夫な気がします。
(追記:最大30x30を想定しておられるとのことでまともに全探索するアプローチでは不十分と思いました)

---

追記：VBAではありませんが、ご要望があったので全探索をするPythonコードを参考までに挙げておきます。
適当に書いたままなのであまりよいコードではありませんが。

import numpy as np
import functools
import operator

A1=[
    [0,1,1,1,1,0,0,0,0,0,0],
    [1,0,0,0,0,1,0,0,0,0,0],
    [1,0,0,0,0,0,1,1,1,0,0],
    [1,0,0,0,0,0,0,0,0,0,0],
    [1,0,0,0,0,0,0,0,0,0,0],
    [0,1,0,0,0,0,0,0,0,2,1],
    [0,0,1,0,0,0,0,0,0,0,0],
    [0,0,1,0,0,0,0,0,0,0,0],
    [0,0,1,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,2,0,0,0,0,0],
    [0,0,0,0,0,1,0,0,0,0,0]
]
A2 = [
    [0,1,0,0,0],
    [1,0,1,0,0],
    [0,1,0,2,1],
    [0,0,2,0,0],
    [0,0,1,0,0]
]

A1 = np.array(A1, dtype=np.uint8)
A2 = np.array(A2, dtype=np.uint8)

def genp():
    # Pを次々に生成するジェネレーター
    # Pは5 x 11とする。1の場所は11個の中から5個選ぶ
    # itertools.permutationsを使えば単純にできるがあえて原始的に実装
    # 組み合わせはTC1 = 11 * 10 * 9 * 8 * 7
    S1 = [i for i in range(11, 6, -1)]
    TC1 = functools.reduce(operator.mul, S1)

    for tc1 in range(TC1):
        i1 = iii(S1, tc1)
        a = np.zeros((5, 11), dtype=np.uint8)
        for i in range(5):
            # print(f'i2[i]={i2[i]}, i1[i]={i1[i]}')
            a[i, i1[i]] = 1
        yield a

def iii(ss, t):
    # t番目の組み合わせのインデックスを生成
    i = [i for i in range(ss[0])]
    r = []
    for s in ss:
        j = t % s
        # print(f'j={j} i={i}')
        r.append(i[j])
        del i[j]
        t //= s
    return r

def main():
    count = 1
    g = genp()
    for p in g:
        A3 = p.dot(A1).dot(p.T)
        if np.all(A2 == A3):
            print(f'number of trials = {count}')
            print('p=')
            print(p)
            break
        count += 1
        # 探索の進捗具合の確認用
        if count % 10000 == 0:
            print(f'count={count}')

if __name__ == '__main__':
    main()

(Python 3.6.5)

Comments

[teratailist]

置換行列のすべての組合せで9!×2 = 725760という数式がありますが、
これはどういった考え方ですか？初歩的な質問のようで申し訳ありません。

実際には、隣接行列A1の大きさは可変で、30行30列になることもあるので、
総当たりだと相当な計算時間が必要になりそうな気がします。

[KSwordOfHaste]

ご質問にあるA1の行列の大きさが11x11で、重みが2の箇所が2つあったため、重み2については交換の組み合わせは2とおりしかなく、他の9個の置換行列の組み合わせは大きさNのときN!だから9!×2と計算しました。
あ・・・A2は5x5だからPは5x11の行列になりますか9!x2にはならないですね。失礼しました。
---
30x30だと素朴に総当たりではやばそうですね。PDFのp.9以降にUllmannの方法以外の最適化（要するに探索空間を狭めるアイデア）があるのでそれらによってどのくらい空間を絞れるかがポイントになると思います。

[KSwordOfHaste]

ちなみにまず総当たりでコードを書いてみてはいかがでしょうか？行列演算が面倒なのでちょっと別の言語で試してみましたが、重みを考慮しない単純な置換行列のしらみつぶし探索でやってみても(自分の論理ではたまたま)4万回弱で発見できました。
---(訂正)--->4万回弱ではなく8万回弱でした。
高度なアルゴリズムにいきなり挑戦するのはなかなか大変なのでまず単純探索論理を書いてみてご質問の例のように小さめのサンプルで解いてみると小手調べになると思います。その上で本格的アルゴリズムへアプローチしてはいかがでしょう。

[teratailist]

ご回答いただきありがとうございます。

参考までに、別の言語で試作したコードを見せていただいてもよろしいでしょうか？
単純な総当たり探索とは言え、まだイメージが持てていないもので…

[KSwordOfHaste]

うーん・・・pythonなのでVBAとは大分違いますし配列演算はnumpyというライブラリーでやってるので仕様をしらないとチンプンカンプンではないかと思います。素朴なコードでしかないので回答に貼るのはかまわないと思いますが・・・

[teratailist]

コードは違えど、ロジックは同じになりますよね。
考え方やpythonのコードを教えていただきたいのですが、いかがでしょうか？
教えていただいたnumpyについては勉強します！

[KSwordOfHaste]

言語が違えば機能も違ってくるのでロジックの考え方も変わってくるものではありますが・・・一応載せておきます。

[teratailist]

ご回答いただきありがとうございます。
解読してみます。

[KSwordOfHaste]

以前のコメントを一つ訂正させていただきます。
元のコードにバグがあって4万回弱で発見したと勘違いしたのですが、本当は８万回弱でした。失礼しました。

[teratailist]

pythonコードを解読する作業を行っております。
main関数の内容をおおよそ理解しました。
今更ながら、置換行列Ｐをどのように生成するかが肝になるようですね。
そこで、genp関数をしっかり理解する必要があると思い、読み込んでいるのですが、
まず、
TC1 = functools.reduce(operator.mul, S1)
では、どのような結果が返されますか？
S1の要素をすべて掛け合わせた数値が一つ返ってくるのかと思ったんですが、
その後に
for tc1 in range(TC1):
が出てきたため、TC1は配列のようだと感じています。
教えてください。よろしくお願いいたします。

[KSwordOfHaste]

TC1はスカラーで11個の中から5個選んだ順列組み合わせ11P5です。

[KSwordOfHaste]

すみません、とんでもない間違いがありました。
上のコメントでようやく気付いたのですが本来の探索空間の大きさは11P5=55440でした。元の論理では11P5x5P5の空間を探索(要するに同じPを複数回生成してしまっていた)ためえらく無駄な探索になってます。
コード例を変更しておきます。変更後の論理では発見に要する試行回数は53791になりました。

[KSwordOfHaste]

ちなみに単純な全探索の空間の大きさはA1のサイズがn, A2のサイズがmならnPmになると思うのでnの大きさはもちろんですがmの大きさに非常に大きな影響を受けると思います。本件のようにm=5ぐらいだとn=30でも探索空間はたかだか10^7程度なので一台のPCでも全探索でOKと思います。しかしmがnに近くなると全世界のPCを投入しても数億年というオーダーになるので相応の最適化が必要と思います。

[teratailist]

ご助言いただきありがとうございます。

想定しているmの最大値は今の所15ですから、
30P15 ≒ 2*(10^20)回の計算が必要になりますね。
30P15回、30行15列のゼロ行列をつくるだけのプログラムを今、VBAで回していますが、
10分間以上実行させても終わっていません。
実際のシステムでは、最大で30回以上この探索を繰り返す(30P15でないにしろ、30P7や30P3など)ので
最適化は必須だと思いました。

[teratailist]

しかし、まずは例にあげている11*5の置換行列を求めるプログラムをVBAで作ろうと思います。

ところで、ご提示いただいたpythonコードのiii関数
t //= s
はどういう処理ですか？

[KSwordOfHaste]

t = t \ s
です。すみませんがこのレベルの質問はご容赦ください。pythonの指導になってしまいます。それは本サイトで本来すべきことではないと思います。

[teratailist]

申し訳ありません。
以降気を付けます。
お陰様でおおよそ理解できました。

Answer 2

お世話になっております。

ご提示いただいたpythonコードを参考に、VBAでも欲しい動作を得ることができました。
ご協力いただいた皆様ありがとうございました。

Option Explicit

Sub main()
    Dim A1  '(1 To 11, 1 To 11)
    Dim A2  '(1 To 5, 1 To 5)
    Dim A3  '(1 To 5, 1 To 5)
    Dim MatT
    Dim MatX
    Dim A1N As Integer
    Dim A2n As Integer
    Dim i, j As Long
    Dim P() As Integer
    Dim prm As Long
    Dim count As Long
    Dim count_ As Long
    Dim i1
    Dim prmttn

    A1 = [{0,1,1,1,1,0,0,0,0,0,0;1,0,0,0,0,1,0,0,0,0,0;1,0,0,0,0,0,1,1,1,0,0;1,0,0,0,0,0,0,0,0,0,0;1,0,0,0,0,0,0,0,0,0,0;0,1,0,0,0,0,0,0,0,2,1;0,0,1,0,0,0,0,0,0,0,0;0,0,1,0,0,0,0,0,0,0,0;0,0,1,0,0,0,0,0,0,0,0;0,0,0,0,0,2,0,0,0,0,0;0,0,0,0,0,1,0,0,0,0,0}]
    A2 = [{0,1,0,0,0;1,0,1,0,0;0,1,0,2,1;0,0,2,0,0;0,0,1,0,0}]
    
    A1N = 11
    A2n = 5
    
    ReDim P(1 To A1N, 1 To A2n)
    
    count = 0
    
    ' 11 * 5のあらゆる置換行列を生成・探索
    prmttn = Permutation(A1N, A2n)
    For prm = 1 To prmttn
        i1 = iii(prm, A1N, A2n)
        
        ' 11 * 5のゼロ行列を生成
        For i = 1 To A1N
            For j = 1 To A2n
                P(i, j) = 0
            Next j
        Next i
        
        ' 置換行列の1の要素を設定
        For i = 1 To A2n
            P(i1(i), i) = 1
        Next i
        
        ' PT * A1 * P (5 * 5)を計算
        MatT = MatrixT(P)
        MatX = MatrixCross(MatT, A1)
        A3 = MatrixCross(MatX, P)
        
        ' A2, A3を比較して、一致する要素をカウント
        count_ = 0
        For i = 1 To A2n
            For j = 1 To A2n
                If A2(i, j) = A3(i)(j) Then
                    count_ = count_ + 1
                End If
            Next j
        Next i
        ' A2, A3の要素がすべて一致すれば、"Pあり"と表示
        If count_ = A2n * A2n Then
            Debug.Print count
            Debug.Print "Pあり"
            
            Exit For
        End If
        
        count = count + 1
        
        ' 探索の進捗状況の確認
        If count Mod 10000 = 0 Then
            Debug.Print count
        End If
        
    Next prm
    
End Sub

' 行列同士の掛け算
Function MatrixCross(ByVal m1, ByVal m2)
    Dim ans
    Dim i, j, k
    Dim sum_
    
    ans = CreateMatrix(UBound(m1), UBound(m2, 2))
    For i = 1 To UBound(ans)
        For j = 1 To UBound(ans(2))
            sum_ = 0
            For k = 1 To UBound(m1(2))
                sum_ = sum_ + m1(i)(k) * m2(k, j)
            Next k
            ans(i)(j) = sum_
        Next j
    Next i
    MatrixCross = ans
    
End Function

' 行列の転置
Function MatrixT(ByVal m)
    Dim ans
    Dim i, j
    
    ans = CreateMatrix(UBound(m, 2), UBound(m))
    For i = 1 To UBound(ans)
        For j = 1 To UBound(ans(2))
            ans(i)(j) = m(j, i)
        Next j
    Next i
    MatrixT = ans
    
End Function

' 任意サイズの行列を作成
Function CreateMatrix(ByVal row_size, ByVal col_size)
    Dim ans, row
    Dim i
    
    ans = Array()
    ReDim ans(1 To row_size)
    For i = 1 To row_size
        row = Array() ' 新しいオブジェクトのインスタンス代入
        ReDim row(1 To col_size)
        ans(i) = row
    Next i
    
    CreateMatrix = ans
    
End Function

' t番目の組合せのインデックスを生成
Function iii(ByVal t As Long, A1N As Integer, A2n As Integer)
    Dim i()
    ReDim i(0 To A1N - 1)
    Dim r()
    Dim k As Integer
    Dim j As Integer
    Dim count As Integer
    count = 0
    
    ' i = 1 ～ A1N
    For k = 0 To A1N - 1
        i(k) = k + 1
    Next k
    
    ReDim r(1)
    For k = A1N To A1N - A2n + 1 Step -1
        count = count + 1
        j = t Mod k
        ReDim Preserve r(count)
        r(count) = i(j)
        Call ArrayRemove(i, j)
        t = t \ k
        
    Next k
    
    iii = r
    
End Function

' 配列の要素の削除
Sub ArrayRemove(ByRef TargetArray, ByVal deleteIndex As Integer)
    Dim i As Integer

    '削除したい要素以降の要素を前につめて上書きコピー
    For i = deleteIndex To UBound(TargetArray) - 1
        TargetArray(i) = TargetArray(i + 1)
    Next i

    '最後の要素を削除
    ReDim Preserve TargetArray(UBound(TargetArray) - 1)

End Sub

' 順列(nPk)の計算
Function Permutation(n, k)
    Dim i
    Dim nPk
    
    nPk = 1
    For i = n - k + 1 To n
        nPk = nPk * i
    Next i
    
    Permutation = nPk
    
End Function

Answer 3

ご提示のコードを動かそうとしましたが、
A1 = [{0,1,1･･･ の処理で既に動きませんが？

このコードはあなたが作ったものですか？
それとも、何かの宿題ですか？

Comments

[teratailist]

申し訳ありません。
テキストファイルでコードを書いて、動作の確認をせずにあげてしまいました。
宿題ではなく、自分で作成したものです。
コードを修正しておきます。

[ExcelVBAer]

コメント欄には修正コード載せなくていいですよ。
コメントが見づらくなっちゃうんで。

あと、↓の判定で、「ElseIf ･･･ = False Then」は、
「Else」でいいのでは？
（２回処理する必要、無いですよね？）
If SubgraphHomotyping(A1, A2) = True Then
Debug.Print "Pあり"
ElseIf SubgraphHomotyping(A1, A2) = False Then
Debug.Print "Pなし"
End If

※「If Pが存在する Then」の箇所も同様

[teratailist]

ご回答いただきありがとうございます。

その通りです。
修正しておきます。