クラスタリング問題

クラスタリング問題とは、データセットをグループに分割する問題です。

グループ内のデータは互いに似ており、異なるグループのデータとは異なる特徴を持っています。

この問題は、機械学習やデータマイニング、ネットワーク解析などの分野でよく用いられます。

クラスタリングによって得られたグループは、同じグループ内のデータが似ているため、そのグループに属するデータの特徴をより深く理解することができます。

また、異なるグループのデータの違いを明確にすることができます。

クラスタリングは、顧客分類、マーケティング分析、医療診断、画像処理、自然言語処理などの分野で幅広く使用されています。

例題

クラスタリング問題の例を１つ挙げます。

Zachary’s karate clubとして知られる、34人のメンバーとそのつながりを示すグラフを用いて、クラスタリングを行います。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

import networkx as nx
import matplotlib.pyplot as plt
from networkx.algorithms import community

# グラフの読み込み
G = nx.karate_club_graph()

# グラフの可視化
pos = nx.spring_layout(G)
nx.draw(G, pos, with_labels=True)
plt.show()

# クラスタリング
communities = community.greedy_modularity_communities(G)

# クラスタリング結果の表示
for i, c in enumerate(communities):
    print("Cluster ", i+1, ": ", c)

# クラスタリング結果をグラフに反映
color_map = []
for node in G:
    for i, c in enumerate(communities):
        if node in c:
            color_map.append(i)
nx.draw(G, pos, node_color=color_map, with_labels=True)
plt.show()

このコードでは、まずグラフを読み込み、可視化します。

次に、community.greedy_modularity_communities()を使用して、グラフをクラスタリングします。

最後に、クラスタリング結果をグラフに反映し、色分けして表示します。

この例では、Zachary’s karate clubグラフは3つのクラスターに分割されます。

各クラスターは、それぞれのメンバーのつながりを考慮して形成されていることがわかります。

[実行結果]

自然災害問題

自然災害に関する問題の一つとして、道路ネットワーク上での交通渋滞の発生があります。

以下のようなシナリオを考えてみましょう：

🔹ある地域において、地震が発生し、複数の道路が通行不能になった。
🔹この地域には複数の車両が存在し、それぞれの車両が目的地に到達するために、どのルートを取るかを決定しなければならない。
🔹通行可能なルートは限られており、道路の通行状況は車両数に応じて変化する。
🔹車両が選択するルートによって、交通渋滞が発生する可能性がある。

このような問題を解くために、以下のような手順でNetworkXを使用することができます。

1. 道路ネットワークのグラフを作成する。
   グラフのノードは交差点や道路の始点・終点とし、エッジは道路とし、エッジの重みはその道路の通行可能度合いとする。
2. 車両の現在地と目的地をノードとして指定する。
3. NetworkXの最短経路アルゴリズムを使用して、車両が最短時間で目的地に到達できる経路を計算する。
4. 計算された最短経路を実際に走行した場合の交通渋滞の発生確率を評価する。

解法・ソースコード

交通渋滞の発生確率を評価する上記の例題をnetworkxを用いて実装します。

具体的には、以下の手順でグラフを作成し、最短経路を計算して、最短経路上での交通渋滞の発生確率を評価しています。

1. グラフを作成する

1
2

G = nx.Graph()
G.add_edges_from([(1, 2), (2, 3), (2, 4), (4, 5)])

まず、NetworkXを用いてグラフを作成します。ここでは、4つのノードと4つのエッジを持つグラフを作成しています。

2. エッジに重みを付ける

1
2

for (u, v) in G.edges():
    G[u][v]['weight'] = random.uniform(0.1, 1.0)

次に、各エッジにランダムな重みを付けます。

ここでは、0.1から1.0の範囲の乱数を用いて重みを設定しています。

3. ノード1からノード5への最短経路を計算する

1
2

path = nx.shortest_path(G, source=1, target=5, weight='weight')
print(f"Shortest path: {path}")

nx.shortest_path関数を用いて、ノード1からノード5までの最短経路を計算します。

weight引数には、エッジの重みを指定しています。

計算された最短経路をpath変数に格納し、print関数を用いて表示しています。

4. 計算された最短経路を実際に走行した場合の交通渋滞の発生確率を評価する

1
2
3
4
5
6
7
8

prob = 1.0
for i in range(len(path) - 1):
    u = path[i]
    v = path[i + 1]
    weight = G[u][v]['weight'] + random.uniform(-0.1, 0.1)
    G[u][v]['weight'] = max(weight - 0.1, 0)  # 通行可能度合いを減少させる
    prob *= weight  # 経路上のエッジの通行可能度合いの積を計算する
print(f"Probability of traffic congestion: {1 - prob:.3f}")

最後に、計算された最短経路上での交通渋滞の発生確率を評価します。

prob変数には、経路上のエッジの通行可能度合いの積を格納します。

forループで、最短経路上の各エッジについて、エッジの通行可能度合いをランダムに減少させ、prob変数にそのエッジの通行可能度合いを乗じます。

最後に、1 - probを計算することで、交通渋滞の発生確率を評価します。

なお、エッジの通行可能度合いの減少幅は、0から0.1の範囲の乱数を用いて決定しています。

以上の手順で、最短経路上での交通渋滞の発生確率を評価することができます。

[実行結果（例）]

Shortest path: [1, 2, 4, 5]
Probability of traffic congestion: 0.953

ノード1からノード5への最短経路は[1, 2, 4, 5]となりました。

そして、この経路上のエッジの通行可能度合いの積を計算することで、交通渋滞の発生確率を評価しました。

結果として、この経路上で交通渋滞が発生する確率は0.953となりました。

最小全域木問題

最小全域木問題 (Minimum Spanning Tree Problem) とは、グラフ理論において、与えられた重みつき無向グラフに対して、そのグラフの全ての頂点を含む木（連結かつ閉路を持たないグラフ）のうち、辺の重みの和が最小になるような木を求める問題です。

問題

以下の図のグラフにおいて、全ての頂点を繋ぐ最小コストの辺集合を求めてください。

     2     3
 A ----- B ----- C
 |   4 / |   5 / |
 |   /   |   /   |
 1 /     4 /     | 1
 | /     | /     |
 |/      |/      |
 D ----- E ----- F
     5     4

解法・ソースコード

最初に、networkxライブラリをnxという名前でインポートします。

次に、空のグラフを作成し、add_weighted_edges_from()メソッドで、頂点間の辺を加えます。

ここで、(‘A’, ‘B’, 2)というように、3つ目の引数で辺の重みを指定しています。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

import networkx as nx

# グラフを作成
G = nx.Graph()
G.add_weighted_edges_from([
    ('A', 'B', 2),
    ('A', 'D', 1),
    ('B', 'C', 3),
    ('B', 'E', 4),
    ('B', 'D', 4),
    ('C', 'E', 5),
    ('C', 'F', 1),
    ('D', 'E', 5),
    ('E', 'F', 4)
])

次に、minimum_spanning_tree()メソッドを使って、最小全域木を求めます。

このメソッドは、引数にグラフを取り、最小全域木を構成する辺集合を返します。

最後に、print()関数を使って、最小全域木を構成する辺集合を表示しています。

edges()メソッドを使うことで、辺集合を取得できます。

1
2
3
4
5

# 最小全域木を求める
mst = nx.minimum_spanning_tree(G)

# 結果を表示
print(mst.edges())

実行すると、以下のように最小全域木を構成する辺集合が表示されます。

[実行結果]

[('A', 'D'), ('A', 'B'), ('B', 'C'), ('B', 'E'), ('E', 'F')]

影響力の分析

影響力の分析についての例題として、Twitter上でのフォロー関係から、影響力の高いユーザーを分析することを考えてみます。

以下のようなツイッターユーザー間のフォロー関係があるとします。

A follows B
A follows C
B follows C
B follows D
C follows A
C follows B
C follows D
D follows A
D follows C

このグラフをNetworkXで表現し、pagerankアルゴリズムを用いて、各ユーザーの影響力を計算するPythonコードを以下に示します。

pagerankアルゴリズムとは、Webページの重要度を計算するアルゴリズムです。

Webページのリンク構造を解析して、そのページが他の多くのページからどの程度の重要度を持つかを計算します。

これを、フォロワー関係からの影響力を分析するのに使ってみます。

[ソースコード]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

import networkx as nx

# ツイッターユーザーのフォロー関係を表現するグラフを作成する
G = nx.DiGraph()
G.add_edges_from([('A', 'B'), ('A', 'C'), ('B', 'C'), ('B', 'D'), ('C', 'A'), ('C', 'B'), ('C', 'D'), ('D', 'A'), ('D', 'C')])

# pagerankアルゴリズムを用いて、各ユーザーの影響力を計算する
pagerank = nx.pagerank(G)

# 影響力が高い順にユーザーを並び替える
sorted_pagerank = sorted(pagerank.items(), key=lambda x: x[1], reverse=True)

# 結果を表示する
for user, score in sorted_pagerank:
    print(f'{user}: {score}')

実行すると、以下のような結果が得られます。

[実行結果]

C: 0.3245609358176832
A: 0.2251463547274389
B: 0.2251463547274389
D: 0.2251463547274389

この結果から、ユーザーCが最も影響力が高いことが分かりました。

問題

以下のグラフにおいて、頂点Aと頂点Eの間にパスが存在するかどうかを判定してください。

  A --- B --- C
  |     |     |
  |     |     |
  D --- E --- F

解法・ソースコード

この問題を解くためには、以下の手順に従います。

1. NetworkXを用いてグラフを作成します。

1
2
3
4
5

import networkx as nx

G = nx.Graph()

G.add_edges_from([('A', 'B'), ('A', 'D'), ('B', 'C'), ('B', 'E'), ('C', 'F'), ('D', 'E'), ('E', 'F')])

2. 頂点Aと頂点Eの間にパスが存在するかどうかを判定します。

NetworkXのhas_path関数を使います。

1
2
3
4

if nx.has_path(G, 'A', 'E'):
    print("頂点Aと頂点Eの間にパスが存在します。")
else:
    print("頂点Aと頂点Eの間にパスは存在しません。")

実行結果は、「頂点Aと頂点Eの間にパスが存在します。」となります。

以上の手順で、グラフにおいて2つの頂点間にパスが存在するかどうかを判定することができました。

最長経路問題

NetworkXを使って簡単な最長路問題を解いてみましょう。

[単純なグラフ]

まず、networkxをインストールして、次のようにインポートします。

1

import networkx as nx

次に、グラフを作成します。ここでは、有向グラフを作成します。

1

G = nx.DiGraph()

グラフにノードを追加します。

1

G.add_nodes_from([1, 2, 3, 4, 5])

ノード間のエッジを追加します。各エッジには、重みがあると仮定します。

1

G.add_weighted_edges_from([(1, 2, 3), (1, 3, 4), (2, 3, 1), (2, 4, 6), (3, 4, 2), (3, 5, 7), (4, 5, 8)])

これでグラフが完成しました。

最長経路問題を解くには、NetworkXのlongest_path関数を使用します。

1
2

longest_path = nx.dag_longest_path(G)
print(longest_path)

これで、最長経路が求められます。

最長経路の出力は、最初のノードから最後のノードまでの最長経路上のノードのリストです。

注意：

この方法は、有向非巡回グラフ（DAG）にのみ適用されます。

循環がある場合は、longest_path関数はエラーを返します。

全ソースコードと結果

上記のコードをまとめると以下のようになります。

[全ソースコード]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

import networkx as nx

# ①グラフの作成
G = nx.DiGraph()

# ②ノードの追加
G.add_nodes_from([1, 2, 3, 4, 5])

# ③エッジの追加
G.add_weighted_edges_from([(1, 2, 3), (1, 3, 4), (2, 3, 1), (2, 4, 6), (3, 4, 2), (3, 5, 7), (4, 5, 8)])

# ④最長経路の計算
longest_path = nx.dag_longest_path(G)

# ⑤結果の出力
print("最長経路: ", longest_path)

実行結果は次の通りです

[実行結果]

最長経路:  [1, 2, 4, 5]

1 ⇒ 2 ⇒ 4 ⇒ 5という最長経路を求めることができました

最短経路問題

NetworkXを使って簡単な最短路問題を解いてみましょう。

例えば、以下のようなグラフがあるとします。

[単純なグラフ]

このグラフにおいて、ノードAからノードDまでの最短経路を求めるソースコードは以下のようになります。

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10
11
12

import networkx as nx

# グラフを定義する
G = nx.Graph()
G.add_edge('A', 'B', weight=5)
G.add_edge('A', 'C', weight=3)
G.add_edge('B', 'D', weight=2)
G.add_edge('C', 'D', weight=3)

# ノードAからノードDまでの最短経路を求める
path = nx.shortest_path(G, 'A', 'D', weight='weight')
print(path)  # 結果：['A', 'C', 'D']

最短経路を求めるために、NetworkXのshortest_path関数を使っています。

またweightパラメータには、辺の重み（コストや所要時間など）を指定しています。

[実行結果]

['A', 'C', 'D']

‘A’ ⇒ ‘C’ ⇒ ‘D’という最短経路を求めることができました

食事プランの最適化

最適化問題の例として食事プランの最適化問題を考えてみます。

この問題では、1週間の食事プランを最適化し、必要な栄養素をバランスよく摂取しながら、予算内で食材を購入することを目的とします。

以下は、PuLPを使用してこの問題を解くためのソースコードです。

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

from pulp import *

# 問題の定義
prob = LpProblem("Meal Planning", LpMinimize)

# 変数の定義
chicken = LpVariable("Chicken", lowBound=0.1)
beef = LpVariable("Beef", lowBound=0.1)
fish = LpVariable("Fish", lowBound=0.1)

# 制約条件の設定
prob += 3 * chicken + 4 * beef + 2 * fish >= 7.5
prob += 2 * chicken + 3 * beef + 1.5 * fish >= 6
prob += 0.5 * chicken + beef + 0.3 * fish >= 1.5
prob += 3 * chicken + 4 * beef + 2 * fish <= 30
prob += 2 * chicken + 3 * beef + 1.5 * fish <= 25
prob += 0.5 * chicken + beef + 0.3 * fish <= 5

# 目的関数の設定
prob += 5 * chicken + 7 * beef + 6 * fish

# 問題の解法
status = prob.solve()

# 結果の出力
print("Status:", LpStatus[status])
print("Optimal Solution:")
for var in prob.variables():
    print(var.name, "=", var.varValue)
print("Total Cost of Ingredients = ", value(prob.objective))

上記のコードでは、変数 chicken、beef、fish はそれぞれ、鶏肉、牛肉、魚介類の購入量を表しています。

prob += 式 不等号 式で制約条件を、prob += 式 で目的関数を設定しています。

制約条件と目的関数の違いは、不等号があるかどうかです。

制約条件は、それぞれ栄養素のバランスを保つための条件です。

例えば、４つ目の制約条件は、タンパク質の摂取量が30グラム以下であることを表しています。

目的関数は、食材の総額を最小化することを目的としています。

最後に、prob.solve() で問題を解き、結果を出力しています。

結果は、各変数(食材)の最適値と、食材の総額です。

このコードを実行すると、以下のような結果が得られます。

[実行結果]

Status: Optimal
Optimal Solution:
Beef = 1.8833333
Chicken = 0.1
Fish = 0.1
Total Cost of Ingredients =  14.2833331

これは、最適な食事プランが、牛肉1.88ポンド、鶏肉0.1ポンド、魚介類0.1ポンドであることを示しています。

また、このプランの食材の総額は14.28ドルであることがわかります。

スケジューリング最適化

スケジューリング最適化の例題を考えます。

複数のジョブを複数のマシンで実行する場合のスケジューリングを最適化する問題です。

各ジョブには処理を実行するために必要な時間があり、各マシンには同時に実行できるジョブの数に制限があります。

目的は、すべてのジョブを完了するのに必要な最小時間を見つけることです。

例えば、4つのジョブと2つのマシンがあり、以下のような情報が与えられたとします。

解き方・ソースコード

この問題をPuLPで解くには、以下のようなソースコードになります。

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

import pulp

# ①問題を定義する
problem = pulp.LpProblem("Job Scheduling", pulp.LpMinimize)

# ②各ジョブをどのマシンで実行するかを決定する変数を作成する
# ジョブAをマシン1で実行する場合、変数A_1の値は1になる
# ジョブAをマシン2で実行する場合、変数A_2の値は1になる
A_1 = pulp.LpVariable("A_1", lowBound=0, cat='Binary')
A_2 = pulp.LpVariable("A_2", lowBound=0, cat='Binary')
B_1 = pulp.LpVariable("B_1", lowBound=0, cat='Binary')
B_2 = pulp.LpVariable("B_2", lowBound=0, cat='Binary')
C_1 = pulp.LpVariable("C_1", lowBound=0, cat='Binary')
C_2 = pulp.LpVariable("C_2", lowBound=0, cat='Binary')
D_1 = pulp.LpVariable("D_1", lowBound=0, cat='Binary')
D_2 = pulp.LpVariable("D_2", lowBound=0, cat='Binary')

# ③目的関数を定義する（すべてのジョブが完了するのに必要な時間を最小化する）
problem += 2*A_1 + 2*A_2 + 3*B_1 + 3*B_2 + 2*C_1 + 2*C_2 + 4*D_1 + 4*D_2

# ④制約条件(1)を定義する（各ジョブは1つのマシンで実行されなければならない）
problem += A_1 + A_2 == 1
problem += B_1 + B_2 == 1
problem += C_1 + C_2 == 1
problem += D_1 + D_2 == 1

# ⑤制約条件(2)を定義する（各マシンは同時に実行できるジョブの数に制限がある）
problem += A_1 + B_1 + C_1 + D_1 <= 2
problem += A_2 + B_2 + C_2 + D_2 <= 2

# ⑥問題を解く
status = problem.solve()

# ⑦結果を表示する
print("Status:", pulp.LpStatus[status])
print("Minimum Time:", pulp.value(problem.objective))
print("A_1:", pulp.value(A_1))
print("A_2:", pulp.value(A_2))
print("B_1:", pulp.value(B_1))
print("B_2:", pulp.value(B_2))
print("C_1:", pulp.value(C_1))
print("C_2:", pulp.value(C_2))
print("D_1:", pulp.value(D_1))
print("D_2:", pulp.value(D_2))

①problemオブジェクトを作成し、最小化問題であることを示します。

②各ジョブをどのマシンで実行するかを決定する変数を作成します。
　各変数は、0または1の2値をとります。

③目的関数を定義します。
　すべてのジョブが完了するのに必要な時間を最小化するように設定しています。

④制約条件(1)を定義します。
　各ジョブは1つのマシンで実行されなければならず、各マシンは同時に実行できます。

⑤制約条件(2)を定義します。
　ジョブの数に制限があることを示します。

⑥問題を解きます。

⑦結果を表示します。
　最適なスケジュールを決定するための変数の値と、最小化された時間を表示します。

このコードを実行すると、以下のような結果が得られます。

[実行結果]

Status: Optimal

Minimum Time: 11.0

A_1: 1.0

A_2: 0.0

B_1: 1.0

B_2: 0.0

C_1: 0.0

C_2: 1.0

D_1: 0.0

D_2: 1.0

この結果から、最小時間は11であり、次のようなスケジュールが最適であることがわかります。

🔹ジョブAをマシン1で、ジョブBをマシン1で、ジョブCをマシン2で、ジョブDをマシン2で実行する。

マーケティング最適化問題

簡単なマーケティング最適化問題を考えます。

この問題では、広告費用を最小限に抑えつつ、与えられた広告予算内で最大の利益を得るために、どの製品にどの程度の広告予算を割り当てるかを決定する必要があります。

例：

ある企業が3つの製品を販売しています。

各製品には、それぞれ異なる利益と広告費用があります。

企業は、合計で10,000ドルの広告予算を持っています。

どの製品にどの程度の広告予算を割り当てるべきかを決定して、利益を最大化する必要があります。

各製品の利益と広告費用は以下のとおりです。

解き方・ソースコード

この問題をPuLPを使って解決するには、以下のようなスソースコードになります。

[Google Colaboratory]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

import pulp

# ①問題を初期化
problem = pulp.LpProblem("マーケティング最適化問題", pulp.LpMaximize)

# ②変数を定義
a = pulp.LpVariable("Aの広告費用", lowBound=0)
b = pulp.LpVariable("Bの広告費用", lowBound=0)
c = pulp.LpVariable("Cの広告費用", lowBound=0)

# ③目的関数を定義
problem += 300*a + 200*b + 150*c, "利益の合計"

# ④制約条件を定義
problem += a + b + c <= 10000, "広告予算の合計"
problem += a <= 2000, "Aの広告費用の上限"
problem += b <= 1500, "Bの広告費用の上限"
problem += c <= 1000, "Cの広告費用の上限"

# ⑤問題を解く
status = problem.solve()

# ⑥結果を出力
print(f"ステータス: {pulp.LpStatus[status]}")
print(f"最大利益: {problem.objective.value()}")
print(f"Aの広告費用: {a.value()}")
print(f"Bの広告費用: {b.value()}")
print(f"Cの広告費用: {c.value()}")

①問題を初期化

pulp.LpProblemを使って問題の名前を指定し、最大化問題であることを指定します。

②変数を定義

pulp.LpVariableを使って各製品の広告費用を表す変数a、b、cを定義します。

また、lowBound引数を使って変数の下限を0に設定します。

③目的関数を定義

各製品の利益を表す定数と広告費用を表す変数をかけた値を足し合わせて、最大化したい利益の合計を定義します。

④制約条件を定義

各制約条件を表す不等式をproblem +=で連結していきます。

制約条件は、広告予算の合計が10000ドル以下、各製品の広告費用が上限に収まることです。

⑤問題解決

problem.solve()を使って問題を解きます。

⑥結果出力

解のステータスや各製品の広告費用、最大利益を出力します。

pulp.LpStatus[status]を使って、解のステータスをテキストで表示します。

また、各変数の値はvalue()を使って取得できます。

最大利益はproblem.objective.value()で取得できます。

[実行結果]

ステータス: Optimal

最大利益: 1050000.0

Aの広告費用: 2000.0

Bの広告費用: 1500.0

Cの広告費用: 1000.0

この結果から、最大利益は1050000ドルであり、A製品には広告費用2000ドル、B製品には広告費用1500ドル、C製品には広告費用1000ドルを割り当てることが最適であることがわかります。

また、制約条件がすべて満たされていることも確認できます。