計測自動制御学会東北支部 第 345 回研究集会(2023.11.07)
資料番号 345-1

HMLasso に対する最適化手法の検討
Examination of optimization methods for HMLasso〇佐々木大夢*，松下慎也*，徐粒*〇Hiromu Sasaki*, Shin-ya Matsushita*， Li Xu**秋田県立大学システム科学技術研究科*Akita Prefectural University
キーワード：HMLasso(Lasso with High Missing rate)，
FISTA (fast iterative shrinkage-thresholding algorithm)，交互方向乗数法(ADMM)
連絡先：〒015-0055 秋田県由利本荘市土谷字海老ノ口 84-4 秋田県立大学 システム科学技術研究科
総合システム工学専攻 制御工学研究室 佐々木大夢 Tel.:(0184)27-2511，
E-Mail:m25n009@akita-pu.ac.jp_____________________________________(1)を解くことはできない．

1. はじめにこの問題に対して先行研究では，(1)に部分𝑋 ∈ ℝ𝑚×𝑛 ，𝑦 ∈ ℝ𝑚 とする．本研究ではラッソ

問題を追加した HMLasso が提案されている

回帰(Lasso)と呼ばれる最適化問題を考える．
1(1)min𝑛 ‖𝑦 − 𝑋𝛽‖22 + 𝜆‖𝛽‖1𝛽∈ℝ 2𝑛
ここで，‖⋅‖2 はℓ2 ノルム，‖⋅‖1はℓ1 ノルム，3)．HMLasso の詳細は 3 節で説明する．先行研究では HMLasso の部分問題を解くために交互方向乗数法(Alternating direction method of

λ > 0とする．適切なλを設定した(1)を解く事

multipliers : ADMM)4)という最適化手法が利用で，成分のほとんどが 0 となる解(スパース解)されていた．一方，ADMM の他にも Fastを得ることができる．(1)は曲線フッティング

Iterative Shrinkage-Thresholding Algorithm

やノイズ処理など広く応用されている 1,2)．

（FISTA）5)と呼ばれる最適化手法が知られてこの問題でカギとなるのが，第一項に含まれいる．FISTA は HMLasso の部分問題に対してる行列𝑋である．実測で得られたデータを𝑋にも適応することができる為，ADMM と FISTA代入するが，測定ミスやセンサの故障などによの有効性を数値実験によって比較する．って，欠損値が多く含まれる場合，そのまま

1

この節では，次節で紹介する HMLasso と最1(2)min𝑛 𝛽 𝑇 𝑆𝛽 + 𝜌𝑇 𝑦 + 𝜆‖𝛽‖1𝛽∈ℝ 211
ここで，𝑆 = 𝑛 𝑋 𝑇 𝑋，𝜌 = 𝑛 𝑋 𝑇 𝑦とする．定義より適化問題を解くアルゴリズムを説明するうえで𝑆は半正定値となる．行列𝑋に欠損が含まれる場必要な数学的な概念を定義する．本節を通して合，式(2)の𝑆，𝜌を，それぞれ以下で定義される

𝑥 ∈ ℝ𝑛，𝐴, 𝐵 ∈ ℝ𝑚×𝑛 ，𝐶 ∈ ℝ𝑛×𝑛 と仮定する．

𝑆 𝑝𝑎𝑖𝑟 ，𝜌𝑝𝑎𝑖𝑟 に変更する．

2. 準備

2.1

𝑆 𝑝𝑎𝑖𝑟 =本研究で使われている数学的な記号

説明
実数全体
n 次元 Euclid 空間ℝ𝑚×𝑛n 行 m 列の行列の集合

𝐴𝑇

行列 A の転置行列していない行の添字集合を示しており，𝑛𝑗𝑘 は集合𝐼𝑗𝑘 の要素の個数を表している．しかし，𝑆 𝑝𝑎𝑖𝑟が半正定値行列とは限らない為，問題(2)の解の存在が保証されないという問題が生じる．

𝒗𝑇 𝐶𝒗 ≥ 0 (∀𝑣 ∈ 𝑅 𝑛 )3. HMLassoℓ1 ノルム𝑛‖⋅‖1‖𝑥‖1 = ∑上記の問題を解決するため，先行研究では|𝑥𝑖 |以下の HMLasso が提案された 3)．

𝑖=1

フロベニウスノルム‖⋅‖𝐹𝐴⊙𝐵(4)

ここで、𝐼𝑗𝑘 は、X の各 j 列と k 列の両方が欠損半正定値行列𝐶𝐶≻01∑ 𝑋𝑖𝑗 𝑦𝑖𝑛𝑗

𝜌𝑝𝑎𝑖𝑟 =𝑖∈𝐼𝑗𝑛ℝ(3)𝑖∈𝐼𝑗𝑘

Table1 本論文で使用する記号
記号ℝ1∑ 𝑋𝑖𝑗 𝑋𝑖𝑘𝑛𝑗𝑘‖𝐴‖𝐹 = √∑𝑗,𝑘1
𝑇min𝑛 𝛽 𝑇 𝜎̂𝛽 + 𝜌𝑝𝑎𝑖𝑟 𝑦 + 𝜆‖𝛽‖1(5)𝛽∈ℝ 22
̂ = arg min‖𝑊 ⊙ (𝜎 − 𝑆 𝑝𝑎𝑖𝑟 )‖𝐹 𝑠. 𝑡. 𝜎 ≻ 0 (6)𝜎2𝑎𝑗𝑘𝜎ここで，𝑊𝑗𝑘 =同じ要素同士の掛け算
𝑎11 𝑏11 ⋯ 𝑎1𝑛 𝑏1𝑛⋮⋱⋮()
𝑎𝑚1 𝑏𝑚1 ⋯ 𝑎𝑚𝑛 𝑏𝑚𝑛𝑛𝑗𝑘𝑛である．HMLasso は欠損率が高いデータに対しても精度の高い解を求めることが出来る．先行研究では部分問題(6)に対して ADMM が適用されていたが，FISTA 等の他の最適化手法を適用した場合の有効性の検証

同じ要素同士の割り算
𝐴⊘𝐵

𝑎11𝑏11( ⋮𝑎𝑚1𝑏𝑚1⋯⋱⋯

𝑎1𝑛𝑏1𝑛が行われていなかった．

⋮ )

𝑎𝑚𝑛
𝑏𝑚𝑛

𝑃𝑆𝐷

半正定値行列の集合

𝑃𝑃𝑆𝐷

PSD の上への射影

4. ADMM と FISTA
この節では，部分問題(6)の解を近似する 2 つの最適化手法を紹介する．
問題(6)に対する ADMM は Algorithm1 で与2.2 欠損値を含んだラッソ回帰問題えられる．

(1)は以下のように表現できる．

2

̂ = 𝑍𝜆̂𝑍 𝑇𝐷(10)Algorithm1 ADMM4)

5. 数値実験

Input : 𝑆 𝑝𝑎𝑖𝑟 , 𝑊, 𝜇, 𝐵1 , 𝐶1for 𝑘 = 1,2, ⋯ 𝒅𝒐
𝐴𝑘+1 = 𝑃𝑃𝑆𝐷 (𝐵𝑘 + 𝑆 𝑝𝑎𝑖𝑟 − 𝜇𝐶𝑘 )𝐵𝑘+1 = (𝐴𝑘+1 − 𝑆 𝑝𝑎𝑖𝑟 + 𝜇𝐶𝑘 )⊘ (𝜇𝑊 ⊙ 𝑊 + 𝐼)
1𝐶𝑘+1 = 𝐶𝑘 + (𝐴𝑘+1 − 𝐵𝑘+1 − 𝑆 𝑝𝑎𝑖𝑟 )この説では，問題(6)の部分問題に対して，交互方向乗数法と FISTA を適用させ結果を比較する．数値実験で用いた計算環境を Table2 に示す．

𝜇end for
Output : 𝐴𝑘

Table2 数値実験で用いた計算環境
OS

FISTA は Algorithm2 で与えられる．

メモリ
CPU

Algorithm2 FISTA5)

Input : 𝑆𝑝𝑎𝑖𝑟

Windows10

使用ソフト

8.00 GBIntel(R) Core(TM) i5-7500
MATLAB

R2023a, 𝑊, 𝜇, 𝐴1 , 𝐵1for 𝑘 = 1,2, ⋯ 𝒅𝒐
𝐴𝑘+1 = 𝑃𝑃𝑆𝐷 (𝐵𝑘 − 𝜇𝑊⨀𝑊⨀(𝐴𝑘 − 𝑆 𝑝𝑎𝑖𝑟 ))t 𝑘+1 = 1 +今回考える問題設定は，以下のとおりである．
𝑋 ∈ ℝ10000×100 とする．ここで𝑋は正規分布に従1+√1+4𝑡𝑘22うように生成する．この行列𝑋をランダムに全体

𝑡 −1の 50％欠損させて，式(4)の𝑆 𝑝𝑎𝑖𝑟 を求める．求𝐵𝑘+1 = 𝐴𝑘+1 + ( 𝑡𝑘 )(𝐴𝑘+1 − 𝐴𝑘 )𝑘+1

end for

めた𝑆 𝑝𝑎𝑖𝑟 を交互方向乗数法(ADMM)と FISTA

Output : 𝐴𝑘を適用し，実行回数ごとに近似解との差分の総和を求めた．これを，各アルゴリズム 10 回初期

Algorithm1 及び Algorithm2 では，𝑃𝑃𝑆𝐷 (半正定点を変え同じ数値実験をし，平均を求めグラフ値の集合(PSD)へ射影)が必要となる．𝐷 ∈ ℝ𝑛×𝑛̂ = 𝑃𝑃𝑆𝐷 (𝐷) とおく.この時， 𝐷̂ は次のに対して, 𝐷にする．初期点は，ランダムな実数としている．

(7)から(10)を計算することで得られる 6)．
(𝐷 + 𝐷 𝑇 )(7)
𝐹=2(8)𝐹 = 𝑍𝜆𝑍 𝑇
式(8)の右辺は行列𝐹の固有値分解を表しており，𝜆は行列𝐹の固有値を対角成分，それ以外に 0 を代入した行列，𝑍は直交行列を表している．また𝜆̂𝑖 (𝑖 = 1,2, … , 𝑛)を次のように設定する．
𝜆𝑖𝜆𝑖 ≥ 0(9)0𝜆𝑖 < 0行列𝜆̂は対角成分が𝜆̂𝑖 ，それ以外に 0 を代入した

Fig.1 部分問題(13)に対して，交互方

𝜆̂𝑖 = {向乗数法(ADMM)，FISTA を実行，解と近似解の距離の平均を比較した行列とし，以下の式(10)の右辺の行列を構成す̂ が得られる．即ち
ると𝐷

結果

3

*http://cvxr.com/

2)

近似解にどれだけ値を近づけることができるか

鈴木譲:スパース推定 100 問 with Python，613,共立出版 (2021)を確認するために，試行回数は 10000 回として3)いる．差分として用いる近似解は，CVX*という

M. Takada, H. Fujisawa and T. Nishikawa:

凸最適化問題を解くソフトウェアを用いて求め

HMLasso: Lasso with High Missing Rate,た値を使用している．結果を Fig.1 で示す.

IJCAI, 3541-3547(2019)4)

Fig.1 から，10000 回 ADMM と FISTA を行

S. Boyd, N. Parikh, E. Chu, B. Peleato, and J.
Eckstein: Distributed optimization andったところ，どちらとも近似解に収束すること

statistical learning via the alternatingが確認できる． 収束の速さは，ADMM のほう

direction method of multipliers. Foundationsが速く近似解に収束することがわかった．それ

and Trends in Machine Learning,1－ぞれの手法を 10000 回実行した平均時間比較し122(2011)

た結果を Table3 に示す．こちらも ADMM の方5)が速いことがわかった.

A. Beck and M. Teboulle: A fast iterative
shrinkage-thresholding algorithm for linear
inverse problems, SIAM J. Imaging Sciences,

Table3 10000 回行った際の実行時間の平均

183-202 (2009)

Algorithm

CPU time[s]

ADMM

28.09

Projection Methods,78-79, Society forFISTA

31.45

Industrial and Applied Mathematics (2011)

6)

6. 結論
本実験では，数値実験でそれぞれのアルゴリズムの有効性を比較することで，交互方向乗数法(ADMM)の方が早く，近似解に収束していることがわかった．
今後の課題は 2 点ある．1 点目は，HMLassoは 2 つの問題を解くことで，解を得ることができるため，計算時間などを考慮するとそれぞれがもっと早く近似解にたどりつく必要がある．
式(6)を解くアルゴリズムはほかにもあるため，他のアルゴリズムを検討し，さらに効率の良い手法を見つけることである．2 つ目は，HMLassoを工場などの実問題に応用させ，有効性を検証することである．

参考文献
1)

永原正章:スパースモデリング―基礎から動的 シ ス テ ム の 応 用 — ， 15 － 71 ， コ ロ ナ 社(2017)

4

R. Escalante and M. Raydan: Alternating