行列計算を用いた機械学習法
今倉暁

観測されたデータや現象を教師データとし，その振る舞いをモデルを用いて再現することで未知のデータの振る舞いを予測する “教師あり機械学習” は，データ解析における重要な問題設定であり，回帰問題や分類問題などに対して，現在さまざまな機械学習アルゴリズムが提案されている．本稿では，特に行列計算に基づく機械学習法として，Ridge 回帰や Lasso などの二乗誤差の最小化に基づく手法，および解析性能の向上のための技術であるカーネル法や次元削減法の概略について紹介する．また，アルゴリズムの実装の際に行列計算の観点からの注意すべき点について数値実験を交えて紹介する．

キーワード：行列計算，機械学習，効率的実装法，カーネル法，次元削減法

1.

はじめに演繹的アプローチに基づく数値解析は，データや現習法が設計される．各ステップの選択肢としては，たとえば下記のようなものが挙げられる．

Step.1 モデルの選択

象の振る舞いを記述する支配方程式を立て，その求解• 線形／非線形回帰モデル

によりデータや現象の原理を厳密に解明することを目•（ディープ）ニューラルネットワーク的とする．一方データ駆動型アプローチに基づくデー• 決定木

タ解析は，データや現象の振る舞いの結果を近似的に• ランダムフォレスト再現することを目的とし，近年活発に研究されさまざ

Step.2 目的関数の定義

まな分野で利活用されている．特に，観測されたデー• 二乗誤差タや現象を “教師データ” とし，その振る舞いをモデル• 相対エントロピーを用いて再現することで，未知のデータの振る舞いを• 交差エントロピー

予測する “教師あり機械学習” は，回帰問題や分類問題などの問題設定に対して，重要な役割として現在さまざまなアルゴリズムが提案されている．
教師データを xi ，対応する正解値を yi とする．教師• 正則化項
Step.3 モデルの最適化•（確率的）勾配降下法•（準）ニュートン法

あり機械学習は，教師データセット {xi , yi }i=1,2,...,n• 行列分解：QR 分解，特異値分解を再現するモデル yi ≈ f (xi , w) (i = 1, 2, . . . , n) を

• Krylov 部分空間法構築することで，テストデータ xtest に対する正解値

本稿では，特に行列計算に基づく機械学習法として，を予測することを目的とする．ここで，f はモデル関Ridge 回帰や Lasso などの二乗誤差の最小化に基づく数，w は重みパラメータ，n は教師データ数である．

手法，および解析性能の向上のための技術であるカー

具体的には以下の四つのステップからなる．

ネル法や次元削減法の概略について紹介する．また，

Step.1 モデル f (xi , w) の選択・設計

アルゴリズムの実装の際に行列計算の観点からの注意

Step.2 目的関数 D({f (xi , w), yi }i=1,2,...,n ) の定義すべき点について数値実験を交えて紹介する．

Step.3 モデルの最適化：
wopt = arg minw D({f (xi , w), yi }i=1,2,...,n )Step.4 テストデータの解析：ytest = f (xtest , wopt )モデル f の選択，目的関数の定義，モデルの最適化法を具体的に設定することで，各種の教師あり機械学
いまくら あきら
筑波大学システム情報系〒 305–8573 茨城県つくば市天王台 1–1–1

imakura@cs.tsukuba.ac.jp

310（14）2. 最小二乗法に基づく解法
データの特徴量数を m，データサンプル数を n とし，教師データセットをX = [x1 , x2 , . . . , xn ] ∈ Rn×m ,y = [y1 , y2 , . . . , yn ] ∈ Rnと置く．以下では，最小二乗法に基づく解法の概要について記す．
オペレーションズ・リサーチ

2.1 最小二乗法
最小二乗法は最も単純な機械学習法の一つであり，モデル関数 f として線形回帰モデル

y ≈ f (X, w) = Xw,w = [w1 , w2 , . . . , wm ] ∈ Rmを用い，誤差として二乗誤差

D(f (X, w), y) = y − Xw22 =2(yi − xi w)iを用いる．
このとき，モデルの最適化は線形最小二乗問題

min y − Xw22w∈Rm(1)として定式化され，行列 X のサイズやランクに注意する必要はあるが，行列 X の QR 分解や特異値分解を用いて厳密に求解できる．
一方データ解析においては，モデルの最適化は必ずしも厳密に行う必要はなく，近似解で十分である場合が多い．この場合，最小二乗問題 (1) は，Krylov 部分

図1(a) L1 正則化と (b) L2 正則化のイメージ空間法などの反復法や，近似 QR 分解や近似特異値分解などの行列の近似分解を用いて効率的に求解するこである Ridge 回帰 [3], Lasso [4] および Elastic Net [5]とができる．最小二乗問題の数値解法については，たの概要について示す．

とえば文献 [1, 2] を参照されたい．

2.2 正則化

2.3.1 Ridge 回帰
線形最小二乗問題 (1) に L2 正則化を加えた方法を

データ解析では，有限個の教師データを用いてモデルの最適化を行う．このとき，単に設定した誤差の最小化を行うと教師データに過剰適合し，未知のデータに対して性能が劣化する．このような過学習を抑制す

Ridge 回帰と呼ぶ（Tikhonov の正則化とも呼ばれる）．
Ridge 回帰におけるモデルの最適化はmin y − Xw22 + λw22(2)w∈Rmるために，誤差項に加えてモデルの複雑さに対するペの求解により行う．ここで，λ > 0 は正則化の重みパナルティを課すことが一般的である．このペナルティラメータである．

を正則化と呼ぶ．代表的な正則化は以下の通りである．

• L2 正則化：w22 =2i wi• L1 正則化：w1 = i |wi |• L0 正則化：w0 = w の非零要素数• 組み合わせ：（例）L1 正則化と L2 正則化

ベ ク ト ル の 2 ノ ル ム の 性 質 a22 + b22=[a , b ] 22 から，最適化問題 (2) の目的関数は⎡⎤ 2⎤ ⎡ y
X22⎦ w⎦−⎣ √y − Xw2 + λw2 = ⎣ 0λIm2

図 1 に 2 次元データに対する L1 正則化と L2 正則のように変形できる．ここで，Im は m 次の単位行列化のイメージを示す．L2 正則化を加えることで，教師である．このため，最適化問題 (2) は，線形最小二乗データへの過剰適合により重み w1 , w2 が過剰に大き問題 (1) と同様に，行列 [X  ,な値を取ることが抑制される．また，L1 正則化を加え特異値分解などの厳密解法や Krylov 部分空間法などることで，重みにゼロ要素が多くなる（図 1(a) の例での近似解法により求解できる．

．
は w1 = 0 となっている）

2.3 具体的な機械学習法
最小二乗問題 (1) に前述の正則化を加えることで，具体的な機械学習法が設計される．以下では代表的手法

2020 年 6 月号

√λIm ] の QR 分解や

2.3.2 Lasso
線形最小二乗問題 (1) に L1 正則化を加えた方法をLasso (Least Absolute Shrinkage and Selection Operator) と呼ぶ．Lasso におけるモデルの最適化は（15）311min y − Xw22 + λw1(3)w∈Rmの求解により行う．

x = φ(x) = [1, x1 , x2 , x1 x2 , x21 , x22 ] ∈ R6 .
簡単のため正則化項を無視すると，非線形変換され最適化問題 (3) は，直接厳密解を求めることができない．このため，反復法で近似的に求解される．代表的なに対する線たデータ X = [x1 , x2 , . . . , xn ] ∈ Rn×m

形回帰は，反復法として，座標降下法（Coordinate Descent, CDmin y − X w22m
w∈R法）や交互方向乗数法（Alternating Direction Method(4)of Multipliers, ADMM 法）[6] が知られている． MAT-のように定式化され，非線形変換されたデータ X にLAB では ADMM 法が標準解法として用いられている．

が得られる．このような非線形対する重み w ∈ Rm

2.3.3 Elastic Net

変換により判別・分類性能は向上するものの，一般に線形最小二乗問題 (1) に L1 正則化と L2 正則化のm  m であるため，最小化問題 (4) の計算コストは両方を加えた方法を Elastic Net と呼ぶ．Ridge 回帰大きく増大する．と同様に，L2 正則化項については二乗誤差と結合し，一方，重み w をLasso と同様に CD 法や ADMM 法などの反復法によw = X w =り求解する．

wi xi ,iw = [w1 , w2 , . . . , wn ] ∈ Rn

2.4 正解データが多次元の場合
正解データが  (≥ 2) 次元である場合，つまり

y i ∈ R ，について考える．このとき，のように xi の線形結合とすると，最小化問題 (4) はmin y − K w22 ,Y = [y 1 , y 2 , . . . , y n ] ∈ Rn×nw∈RK = X X  ∈ Rn×n

(5)と置くと，線形回帰モデルの重みは m ×  次の行列 Wと変形される．ここで，K = X X  をグラム行列と2ij aij )呼ぶ．高次元データ xi = φ(xi ) を計算することなく，となり，行列の Frobenius ノルム (AF =を用いて，最小二乗問題 (1) はminW ∈Rm×

直接グラム行列 K を計算することができれば，最小化問題 (5) を用いることで，計算の効率化が可能である．

Y − XW 2F

このように元データの非線形変換を行うことなく，と書き換えられる．このとき，正則化項としては，L2 正則化の拡張である W 2F や L1 正則化の拡張として，

W 2,1 =

ij2wijkij = φ(xi ) φ(xj ) = k(xi , xj )形回帰を行う手法をカーネルトリックと呼ぶ．ここで，関数 k(xi , xj ) をカーネル関数と呼ぶ．カーネルトリックを用いることで，計算量の削減の他，無限次元の非解析性能の向上のための技術線形関数を利用することも可能となる．

本節では，解析性能向上のための技術として，カーネル法および次元削減法について記す．
解析性能の改善のため，入力データ xi ∈ R

代表的なカーネルとして，多項式カーネル（c ≥ 0, d：自然数）

3.1 カーネル法

,xi = φ(xi ) ∈ Rm

K = [kij ],のように直接計算し，最小化問題 (5) を解くことで線などが用いられる．

3.

グラム行列 K をmをm>mのように非線形変換することを考える．このような非線形変換は，線形変換のみでは判別・分類できない問題に対して有効である．さまざまな非線形関数 φ(x) が考えられるが，たとえば 2 次元データ x = [x1 , x2 ]

dk(xi , xj ) = x,i xj + c
シグモイドカーネルk(xi , xj ) =1,1 + exp(−γxi xj )ガウスカーネル（RBF カーネル）xi − xj 22k(xi , xj ) = exp −,σ2(6)に対する 2 次の多項式を用いると，以下のようになる．などが知られている．

312（16）

オペレーションズ・リサーチ

3.2 特徴量選択・次元削減法

一方で，上記の方法は m 次元に次元削減する際に機械学習において，数学的には，特徴量の次元数 mm + 1 番目以降の固有ベクトルを捨てるため，有用なが大きいほど教師データにおける最適化問題の関数値情報を捨てている可能性がある．これに対して近年，を小さくすることができる．しかしながら，実用上は複素モーメント型並列固有値解法 [11, 12] のアイディ

特徴量の次元数が大きすぎると過学習などにより解析

アに基づき，m 次元データの構築に際し m + 1 本目性能が低下することが知られており，また解析手法の以降の固有ベクトルも併せて利用する複素モーメント計算コストが増大するといったデメリットがある．

型次元削減法 (Complex Moment-based Supervised

このような観点から，特徴量の次元数を削減し解析手法の計算時間の削減および解析性能の向上を行う手法として，特徴量選択法および次元削減法がある．

3.2.1 特徴量選択法

Eigenmap, CMSE)[13] が提案され，高い認識性能を実現することが示されている．

4. 行列計算の観点での計算の工夫

元データの一部の特徴量を選択的に抽出し低次元デー

近年の計算機においては， CPU における浮動小数タを作成する方法を特徴量選択と呼ぶ．特徴量選択法

点演算の演算速度に対して， CPU とメモリの間でのは，各特徴量の統計量を独立に評価し選択する Filter

データ転送速度が非常に遅い．このため，計算時間は法，選択した特徴量に対する解析モデルの性能により演算回数だけでなくデータ転送量やデータ転送回数に適切な特徴量の組み合わせを選択する Wrapper 法，解強く依存し，使用するデータ量に対して演算回数が多析モデルで計算した重みをもとに特徴量の重要度を算い計算ほど計算性能（単位時間あたりの浮動小数点演出する Embedded 法に分けられる [7]．

算回数）が高くなることが知られている．基本線形代2.3.2 節で示した Lasso は，L1 正則化により重みに数ライブラリ (Basic Linear Algebra Subprograms,ゼロ要素が多くなり，これにより重要度の高い特徴量BLAS) のレベル 1, 2, 3 に対応する代表的な行列計算を選択する Embedded 型の特徴量選択手法としてもである，ベクトルの加減算やノルム計算などのベクト用いられる．

ル演算（行列の加減算も含む），行列ベクトル積，行列3.2.2 行列トレースの最適化に基づく次元削減法
を用い，高次元データ xi ∈ Rm
射影行列 B ∈ Rm×m

を低次元（m < m 次元）に射影する手法

xi = B  xi ,

B ∈ Rm×m

を次元削減法と呼び，教師なし次元削減法である主成分行列積の計算性能は以下の通りである．
行列行列積  行列ベクトル積  ベクトル演算．
この性質に着目し，本稿では，ガウスカーネルのグ
ラム行列 K の計算 (6) および LPP や LFDA などトレースの最適化に基づく次元削減法で用いられる行列分析 (Principal Component Analysis, PCA)，局所性

wij (xi − xj )(xi − xj ) ∈ Rm×m ,保存射影 (Locality Preserving Projections, LPP)[8]S=および教師あり次元削減法であるフィッシャー判別分析

wij = wji

(7)

ij

(Fisher Discriminant Analysis, FDA)[9]，局所フィッシャー判別分析 (Local FDA, LFDA)[10] などが代表の計算について，実装法による計算時間の違いを評価的な手法として知られている．

する．本稿では，例として MATLAB コードを示す．

これらの手法は，各解法においてそれぞれ設定される対称行列 A1 ∈ R

m×m

4.1 グラム行列 K の計算法

および正定値対称行列サンプル数 n，パラメータ σ および教師データ X =A2 ∈ Rm×m を用い，行列トレースの最適化問題
min Tr B  A1 B or max Tr B  A1 Bmm
B∈Rm×B∈Rm×
s.t. Tr B  A2 B = 1[x1 , x2 , . . . , xn ] ∈ Rn×m が n, sigma および 2 次に帰着され，実対称一般化固有値問題

A1 u = λA2 u,λ ∈ R,u ∈ Rm \ {0}の最小／最大から m 個の固有値に対応する固有ベクトルとして求解される．

2020 年 6 月号

元配列 X として与えられたとし，グラム行列 K を保持する 2 次元配列 K を計算する．

4.1.1 ナイーブな実装
ガウスカーネルのグラム行列 K は，式 (6) に基づき下記のように計算される．
K = zeros(n);for j = 1:n
for i = 1:n（17）313

上記の実装法の計算コストの主要部は XX  の行列

K(i,j) =...
exp(-norm(X(i,:)-X(j,:))^2 / sigma^2);
endend

積計算であり，計算量はベクトル演算に基づく実装法とおおむね変わらないものの，計算効率が高く計算時間が短くなることが予想される．

2

上記のナイーブな実装法では，n 回のノルム計算が4.2 行列 S の計算法
特徴量次元数 m，サンプル数 n，重み行列 W = [wij ]計算コストの主要部である．

4.1.2 対称性を利用した実装

および教師データ X = [x1 , x2 , . . . , xn ] ∈ Rn×m がグラム行列の対称性 K = K  を利用すると，グラm, n および 2 次元配列 W, X として与えられたとし，
式 (7) の行列 S を保持する 2 次元配列 S を計算する．

ム行列 K の計算は下記のように書き換えられる．

4.2.1 ナイーブな実装

K = zeros(n);for j = 1:n
K(j,j) = 1;
for i = 1:j-1
K(i,j) =...
exp(-norm(X(i,:)-X(j,:))^2 / sigma^2);
K(j,i) = k(i,j);
end
end

行列 S は，式 (7) に基づき下記のように計算される．
S = zeros(m);for j = 1:n
for i = 1:nxd = X(i,:) - X(j,:);S = S + W(i,j) * xd’ * xd;
end
end

対称性を利用した実装法では，n(n − 1)/2 回のノル計算コストの主要部である．

ム計算が計算コストの主要部である．

4.2.2 対称性を利用した実装

4.1.3 行列行列積に基づく効率的な実装
計算性能の低いベクトル演算に基づく上記の実装法
に対して，計算効率の高い行列行列積に基づく実装を考える．グラム行列 K の計算式 (6) はK = [kij ],kij = exp −gij ,σ2

上記のナイーブな実装法では，n2 回の行列の加算がgij = xi − xj 22のように書くことができる．ここで，gij = xi − xj 22 = xi xi + xj xj − 2xi xj

重みの対称性 wij = wji を利用すると，行列 S の計算は下記のように書き換えられる．
S = zeros(m);for j = 1:n
for i = 1:j-1
xd = X(i,:) - X(j,:);S = S + W(i,j) * xd’ * xd;
endendS = 2 * S;であるため，行列 G = [gij ] は，

g = [x1 22 , x2 22 , . . . , xn 22 ] ∈ Rn(8)および 1 = [1, 1, . . . , 1] ∈ Rn を用い，対称性を利用した実装法では，n(n − 1)/2 回の行列の加算が計算コストの主要部である．行列 S の対称性も利用することでさらなる計算量の削減が可能である．

G = [gij ] = g1 + 1g  − 2XX 

4.2.3 行列行列積に基づく効率的な実装
計算性能の低いベクトル演算に基づく上記の実装法のように計算できる．したがって，ガウスカーネルのに対して，計算効率の高い行列行列積に基づく実装を

グラム行列 K は行列行列積に基づき下記の通り計算さ考える．行列 S の定義式 (7) はれる．

g = sum(X.^2,2);
G = repmat(g,1,n) + repmat(g’,n,1) - 2*X*X’;
K = exp(-G/sigma^2);wij (xi − xj )(xi − xj )
ij 
=2
wij xi x
wij xi x
i −2j

ここで，g = sum(X.^2,2); はベクトル g の計算 (8)= 2X  (D − W )X,

S=iの MATLAB コマンドであり，repmat(g,1,n) および repmat(g’,n,1) は g1

MATLAB コマンドである．
314（18）および 1g

を計算する

W = [wij ],jD = diag(di ),ij

di =wij(9)jのように書き換えることができる．したがって，行オペレーションズ・リサーチ

表1

グラム行列 K の計算時間（秒）S の計算において，行列行列積に基づく実装法は，ベサンプル数 n

クトル演算に基づく実装法と比較して計算量が 1/100

1000

20004000程度に削減されているが，計算量削減効果以上の大幅

7.78
3.85
0.03

40.04
19.95
0.11

202.4199.35
0.43

な高速化を実現している．これは，計算時間が計算量

実装法
ナイーブな実装
対称性を利用した実装
行列行列積に基づく実装

表2

ナイーブな実装
対称性を利用した実装
行列行列積に基づく実装

クトル演算と比較して高速であることに由来する．

5. おわりに

行列 S の計算時間（秒）実装法

だけではなく計算効率に強く依存し，行列行列積がベ本稿では，特に行列計算に基づく機械学習法として，サンプル数 n1000

2000

4000

Ridge 回帰や Lasso などの二乗誤差の最小化に基づく

15.29
8.31
0.01

72.88
36.11
0.03

300.67150.72
0.10

手法，および解析性能の向上のための技術であるカーネル法や次元削減法の概略について紹介した．また，アルゴリズムの実装の際に行列計算の観点からの注意すべき点について数値実験を交えて紹介した．

列 S は行列行列積に基づき下記の通り計算される．

数値実験で示されたように，アルゴリズムの計算時間は計算量だけでなく計算効率に強く依存し，効率的

D = diag(sum(W));

な実装を用いることで，時として大幅な高速化が可能

S = 2 * X’ * (D-W) * X;である．このため，計算効率の高い行列行列積計算を基盤とした実装を行うことが非常に重要である．また，

ここで，D = diag(sum(W)); は対角行列 D の計算(9) の MATLAB コマンドである．

実装を行うことができない場合もあるため，アルゴリ

上記の実装法では，計算量が 1/m 程度に削減されている．さらに，計算コストの主要部は X (D − W )Xの行列積計算であり，計算効率が高く計算時間が短くなることが予想される．

4.3 各実装法の計算時間の評価
ガウスカーネルのグラム行列 K の計算 (6) および行列 S の計算 (7) に対して，実装法による計算時間の違いRを評価する．実験環境は，Windows 10 Pro, Intel

CoreTM i7-10710U CPU @ 1.10GHz, 16GB RAMであり，MATALB2019b を使用した．
グラム行列 K の計算 (6) に対しては，特徴量数をm = 1000 とし，サンプル数を n = 1000, 2000, 4000，σ = 1 とし，データ X は乱数により生成した．一方，
行列 S の計算 (7) に対しては，特徴量数を m = 100とし，サンプル数を n = 1000, 2000, 4000 とし，データ X および重み W は乱数により生成した．
実装法毎のグラム行列 K および行列 S の計算時間をそれぞれ表 1 および表 2 に示す．実験結果から，対称性を利用した実装はナイーブな実装と比較して計算時間が半減していることがわかる．これは，対称性を利用することで計算量が半減していることに由来する．
一方で，行列行列積に基づく実装法は，ベクトル演算に基づく実装法と比較して，特に大きな n に対して大幅に計算時間が減少していることがわかる．なお行列2020 年 6 月号

アルゴリズムによっては行列行列積に基づく効率的なズム選択の段階で計算効率を考慮する必要がある．
各種行列計算のアルゴリズムや高性能実装について
は文献 [1, 2] を参照されたい．
参考文献[1] 日本応用数理学会監修，櫻井鉄也，松尾宇泰，片桐孝洋編，『数値線形代数の数理と HPC』共立出版，2018.
[2] G. H. Golub and C. F. Van Loan, Matrix Computations, Fourth Edition, JHU Press, 2013.
[3] A. E. Hoerl and R. W. Kennard, “Ridge regression:biased estimatin for nonorthogonal problems,” Technometrics, 12, pp. 55–67, 1970.
[4] R. Tibshirani, “Regression shrinkage and selection
via the Lasso,” Journal of the Royal Statistical Society (Series B), 58, pp. 267–288, 1996.
[5] H. Zou and T. Hastie, “Regularization and variableselection via the elastic net,” Journal of the Royal Statistical Society (Series B), 67, pp. 301–320, 2005.
[6] S. Boyd, “Distributed optimization and statisticallearning via the alternating direction method of multipliers,” Foundations and Trends in Machine Learning,
3, pp. 1–122, 2010.
[7] G. Chandrashekar and F. Sahin, “A survey on feature selection methods,” Computers and Electrical Engineering, 40, pp. 16–28, 2014.
[8] X. He and P. Niyogi, “Locality preserving projections,” In Advances in neural information processingsystems, pp. 153–160, 2004.
[9] R. A. Fisher, “The use of multiple measurements intaxonomic problems,” Annals of Human Genetics, 7,
pp. 179–188, 1936.
[10] M. Sugiyama, “Dimensionality reduction of mul-（19）315

timodal labeled data by local Fisher discriminantanalysis,” Journal of Machine Learning Research, 8,
pp. 1027–1061, 2007.
[11] T. Sakurai and H. Sugiura, “A projection methodfor generalized eigenvalue problems using numericalintegration,” Journal of Computational and AppliedMathematics, 159, pp. 119–128, 2003.
[12] A. Imakura, L. Du and T. Sakurai, “Relationships316（20）among contour integral-based methods for solving generalized eigenvalue problems,” Japan Journal of Industrial and Applied Mathematics, 33, pp. 721–750, 2017.
[13] A. Imakura, M. Matsuda, X. Ye and T. Sakurai,“Complex moment-based supervised eigenmap for dimensionality reduction,” In Proceedings of the AAAI
Conference on Artificial Intelligence, 33, pp. 3910–3918, 2019.

オペレーションズ・リサーチ