制御理論的スプラインのスパース表現○永原正章 （京都大学） 山本裕 （京都大学）

Sparse Representation for Control Theoretic Splines∗M. Nagahara (Kyoto University) and Y. Yamamoto (Kyoto University)
Abstract– In this article, we consider sparse representation for control signals which is transmitted througha rate-limited communication channel. The problem is formulated as 1 -penalized control theoretic splines.
Then sparse solution is obtained by 1 -2 optimization which can be eﬀectively solved by iterative shrinkage
method. Design examples are shown to illustrate the eﬀectiveness of the proposed method.
Key Words: Sparse representation, control theoretic splines, 1 -2 optimization

1はじめにで定義する．また行列 Φ に対して，そのノルムを圧縮センシングまたは圧縮サンプリングは，元のデータがスパースであるという仮定を置くことによって，非常に少ない観測データから元のデータをかなり正確に復元できる手法として，近年，とくに信号処理や情報理論の分野で盛んに研究されている 2, 6) ．この手法の特徴は，スパースネスの指標であるベクトルの非ゼロ要素の数（しばしば “0 ノルム” と呼ばれる）の最小化を，その凸緩和である 1 ノルムの最小化に置き換えて逆問題を解く点にある．これは，最小二乗法の正則化項に 1 ノルムを用いる LASSO (Least AbsoluteShrinkage and Selection Operator) と呼ばれるものとも関係がある 6) ．この最適解は，内点法を用いて計算することもできるが，近年，繰り返し縮小法 (iterativeshrinkage) と呼ばれる手法が提案され，より高速に最適解を見つけることが可能となった 1, 8) ．
本論文では，この手法にもとづき，制御入力信号のスパース表現を試みる．リモートコントロールでは，システムへの制御入力信号が帯域制限のある通信路を通ってシステムに加えられる．このとき，この制御入力信号のスパース表現（すなわち，送信信号の多くが 0 であるような表現）が可能であれば，例えばエントロピー
符号化 4) の手法を用いて効率的な情報圧縮が可能となる．この性質にもとづき，システムの出力を 2 ノルムで最適化する制御入力の設計問題に設計パラメータの1 ペナルティを含んだ枠組みでの最適化問題を定式化し，設計パラメータのスパース表現を得ることを提案する．具体的には，S. Sun ら 5) によって提案された制御理論的スプラインの設計において，パラメータに対する 2 ペナルティのかわりに 1 ペナルティを導入し，パラメータのスパース表現を行う．この問題は，1 /2最適化として定式化され，上記で述べた繰り返し縮小法により効率的に最適解を得ることができることを示す．また，設計例により，本手法の有効性を示す．
記法ベクトル v = [v1 , . . . , vn ] ∈ Rn に対して，その 1ノルムおよび 2 ノルムを

v1 :=ni=1|vi |,v2 := n
i=1Φ := max

v2 =0で表す．ここで σmax (Φ) は行列 Φ の最大特異値である．実数 x ∈ R に対して

1, if x ≥ 0,sgn(x) :=−1, if x < 0,(x)+ := max{x, 0}と定義する．区間 [0, T ] (T > 0) 上の 実数値 Lebesgue
2 乗可積分関数の全体を L2 [0, T ] と書く．f, g ∈ L2 [0, T ]に対して，内積を T
f, g :=f (t)g(t)dt0で表す．

2問題設定
次の線形システム P を考える:
ẋ(t) = Ax(t) + Bu(t),P :y(t) = Cx(t),t ∈ [0, ∞).

ここで A ∈ Rn×n , B ∈ Rm×1 and C ∈ R1×n とし，x(0) = 0 とする．また，(A, B) は可制御，(C, A) は可観測とする．このシステムの出力信号 y に対して，次のデータ点を与える:

D := {(t1 , α1 ), (t2 , α1 ), . . . (tN , αN )}.
ここで，各 ti は時刻を表すデータで

0 = t1 < t2 < · · · < tN =: Tとし，その各時刻において，システムの出力 y(ti ) がαi に近くなるように，すなわち y(ti ) ≈ αi となるように制御入力 u を設計したい．
この設計問題に対して，次の 2 ノルム評価関数を考える:

1/2vi2Φv2= σmax (Φ)
v2J2 (u) := λ T02u(t) dt +N
i=1wi (y(ti ) − αi )2 .

(1)

D(λW−1θ

−1+ G0 )gi (t)

uPyFig. 1: Remote control system optimized with J2 (u)in (1). The vector θ is transmitted through a communication channel.
ここで，wi > 0 は各データ点へのフィッティングに対する重みである．この評価関数 J2 (u) を最小化する最適制御入力は，次式で表されることが知られており，制御理論的スプラインと呼ばれる 5, 3) :u(t) =N

θi gi (t).

(2)CeA(ti −t) B,gi (t) :=0,if ti > t,if ti ≤ t(3)で定義され，係数 θ := [θ1 , . . . , θN ] は

θ = (λW −1 + G0 )−1 α(4)で与えられる．ここで，α := [α1 , . . . , αN ] であり，
また⎤⎡
gt1 , gt1  . . . gtN , gt1 ⎥
⎢
..
..
..
G0 := ⎣⎦,.
.
.
gt1 , gtN  . . . gtN , gtN 

W := diag {w1 , . . . , wN }である．リモートコントロールの枠組みで考えると，この制御では，まずデータ D が与えられたとき，(4) によって係数ベクトル θ が計算され，この θ が通信路に送信される．受信側では，(2) により制御入力信号 uが計算される．この状況では，通信路の容量に制約があるのが通常で，θ のデータ量は少ない方が望ましい．
したがって，圧縮センシングの考え方 2, 6) に従い，送信信号 θ のスパース表現を考える．まず，制御入力 u
は，L2 [0, T ] の部分空間
M
VM := u ∈ L2 [0, T ] : u =θi ψi , θi ∈ R . (5)i=1に属する信号と仮定する（M ≥ N とする）．ここで，{ψi } は L2 [0, T ] の線形独立なベクトルとする．このとき，評価関数 J2 (u) の第二項は次のように表現できる:

wi (y(ti ) − αi )2i=1= θ G W Gθ − 2θ G W α + α W α
2= Φθ − v2 .
⎡g1 , ψ1 ⎢..
G := ⎣.
g1 , ψM (6)

上で述べたように，第一項の 1 ノルムはベクトル θ のスパースネスを表す指標である．この評価関数を最小化する最適解を θ∗ とすると，最適制御入力はu∗ (t) =Nθi∗ ψi (t),t ∈ [0, T ]で与えられる．

3関数 gi (t) は

ただし，2J1 (θ) := λθ1 + Φθ − v2 .

i=1i=1Nとし，Φ := W 1/2 G, v := W 1/2 α とする．これより，
係数ベクトル θ に 1 のペナルティを与えた次の評価関数を導入する:

...
..
.
...

⎤gN , ψ1 ⎥..
⎦.
gN , ψM 

1 -2 最適化によるスパース表現(6) の評価関数は θ に関して凸であり，最適解は必ず存在する．しかし，その最適解は (4) のように解析的には表すことができず，最適解を求めるためには繰り返し計算が必要である．この繰り返し計算には，内点法などが従来用いられてきたが，近年，繰り返し縮小法 (iterative shrinkage)1, 8) と呼ばれる手法が提案され，これにより高速な数値計算が可能となった．
繰り返し縮小法のアルゴリズムは以下で与えられる
1, 8): 初期値 θ[0] ∈ RN を与え，β[1] = 1, η[1] = θ[0]とする．定数 c をc > Φ = σmax (Φ)となるように選ぶ．次を繰り返す1 :1 
Φ (v − Φη[k]) + η[k]c

θ[k] = Sλ/c
1 + 1 + 4β[k]2β[k + 1] =2β[k] − 1η[k + 1] = θ[k] +(θ[k] − θ[k − 1])β[k + 1]
k = 1, 2, . . .
関数 Sλ/c は η = [η1 , . . . , ηN ] に対して

⎤sgn(η1 )(|η1 | − λ/c)+⎥⎢..
Sλ/c (η) := ⎣⎦.
⎡sgn(ηN )(|ηN | − λ/c)+で定義される非線形関数である．関数 Sλ/c の各要素に対する非線形関数 sgn(η)(|η| − λ/c) の入出力関係をFig. 2 に示す．条件 c > Φ が成り立つような c > 0を選べば，上のアルゴリズムは任意の Φ と v に対して収束し，limk→∞ θ[k] は (6) を最小化する最適解に一致する 1) ．さらに，最適解への収束の早さは最悪評
価で O(1/k 2 ) であることも知られている 1) ．上記のアルゴリズムは，計算機に容易に実装が可能であり，収束も早いことから，リアルタイムでスパースなベクトル θ を得ることが可能である．
1 繰り返し縮小法にも様々な方法が提案されているが 8) ，ここで用いた方法は FISTA (Fast Iterative Shrinkage-ThresholdingAlgorithm) と呼ばれる方法である 1) ．

θy10.5−λ/cη0

λ/c0-0.5-1-1.5Fig. 2: Nonlinear function θ = sgn(η)(|η| − λ/c)+
Original curve and sampled data0123456Fig. 4: Original curve sin t: (blue line), Output y(t)of 1 -2 optimal input with λ = 0.02 (red line), 2
optimal input with λ = 0.1 (blue dots).

10.8
0.60.4

u50.2
04-0.23-0.4
2-0.6-0.8
-1

1012

3

4

560Fig. 3: Original curve sin t and sampled data D.

4-1-2

設計例

⎡ ⎤0
B = ⎣0⎦ ,
112

3

4

56Fig. 5: 1 -2 optimal input u(t) with λ = 0.02 (redline), 2 optimal input u(t) with λ = 0.1 (blue dots).

線形システム P の実現を⎡⎤
0 1 0A = ⎣0 0 1⎦ ,0 0 00

下で与えられる．

C= 10
0 .とする．すなわち，P (s) = 1/s3 である．データ D のサンプル時刻を ti = 0.1+(i−1)×0.5, i = 1, 2, . . . , 11}とし，データ αi を正弦波 αi = sin ti とする（図 3 を参照）．また，L2 [0, T ] の線形独立なベクトル ψi (t) は，上記の A, B, C および {t1 , . . . , t11 } を使って (3) によって定義される gi (t) と同じものとする2 ．また，(1)に対する λ を 0.1，(6) に対する λ を 0.02 とする．これらの設定のもとで，(2) で計算される， (1) の 2 評価関数を最適化する θ を θ2∗ とおき，また (6) の 1 -2
評価関数を最適化する θ を θ1∗ とおくと，それぞれ以2 (A, B) が可制御かつ (C, A) が可観測ならば，(3) で定義される {gi } は線形独立となることが知られている 3) ．

⎤
0.9919⎢ 4.7303 ⎥⎥⎢⎢ 5.4459 ⎥⎥⎢⎢ 3.9878 ⎥⎥⎢⎢ 1.8978 ⎥⎥⎢⎥θ2∗ = ⎢
⎢ 0.1898 ⎥ ,⎢−0.8228⎥⎥⎢⎢−1.2059⎥⎥⎢⎢−1.0451⎥⎥⎢⎣−0.3015⎦
1.0642⎡

⎤
0⎥
⎢
0⎥
⎢
⎢ 2.0665 ⎥⎥⎢⎢ 3.8879 ⎥⎥⎢⎢ 3.6590 ⎥⎥⎢⎥θ1∗ = ⎢
⎢ 1.1288 ⎥⎥⎢0⎥⎢⎢−1.7358⎥⎥⎢⎢−2.6657⎥⎥⎢⎣−0.3231⎦
1.7432⎡

1 -2 最適化によりスパースな解 θ1 が得られることがわかる．これらのパラメータに対するシステムの出力を Fig. 4 に，入力 u(t) を Fig. 5 に示す．また復元誤
差 e(t) = | sin t − y(t)| を Fig. 6 に示す．1 -2 による最適化によるスパース表現にもかかわらず，2 最適化とほぼ同じ性能を示していることがわかる．
さらに，(6) の 1 -2 評価関数の λ を 0.05 にした場error

error

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0012345Fig. 6: Absolute error with θ1∗ (solid) and θ2∗ (dots).

0123456Fig. 8: Absolute error with θ1∗ (solid) and θ2∗ (dots).

5u5

おわりに

本稿では，リモートコントロールの枠組みで，送信信号（最適制御入力を定義するパラメータ θ）のデータ量を抑えるために，1 -2 最適化の枠組みで問題を定式化し，繰り返し縮小法により効率的に最適解を求めることができることを示した．また設計例を示し，性能をそれほど下げずにある程度のスパース表現が可能なこと，またパラメータのスパースネスと制御性能の間にトレードオフの関係があることを示した．

43210

参考文献

-1-2060

123456Fig. 7: 1 -2 optimal input u(t) with λ = 0.05 (redline), 2 optimal input u(t) with λ = 0.1 (blue dots).
合の最適パラメータ θ1∗ は以下のようになる．

⎤0⎥⎢0⎥⎢⎢ 0.3672 ⎥⎥⎢⎢ 3.2113 ⎥⎥⎢⎢ 3.2282 ⎥⎥⎢⎥.
0θ1∗ = ⎢⎥⎢⎥⎢
0⎥⎢⎢−0.6546⎥⎥⎢⎢−2.2377⎥⎥⎢⎦⎣0⎡1.0940定数 λ の値を大きくしたことで，よりスパースな表現が得られていることがわかる．この θ1∗ を用いた入力u(t)， および復元誤差をそれぞれ Fig. 7，Fig. 6 に示
す． 入力 u(t) は λ = 0.02 の場合に比べて小さくなっていることがわかる．これは，1 ノルムのペナルティをより大きくとったことに対応している．さらに，このペナルティを大きくすることにより，復元誤差が上昇していることがわかる．この例題から，入力信号のパラメータ θ のスパースネスと復元誤差（性能）との間にはトレードオフの関係があることがわかる．

1) A. Beck and M. Teboulle, A fast iterativeshrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imaging Sciences, vol. 2,No. 1, 183/202 (2009)2) E. J. Cnades and M. B. Wakin, An introduction
to compressive sampling, IEEE Signal Processing
Magazine, vo. 25, 21/30 (2008)3) M. Egerstedt and C.F. Martin, Control TheoreticSplines, Princeton University Press (2010)4) G. A. Jones and M. Jones, Information and Coding Theory, Springer (2000)
5) S. Sun, M. B. Egerstedt, and C. F. Martin, Control theoretic smoothing splines, IEEE Trans. onAutomatics Control, Vol. 45, No. 12, 2271/2279(2000)
6) 田中，圧縮センシングの数理，IEICE Fundamental
Review, Vol. 4, No. 1, 39/47 (2010)
7) K. L. Teo and X. Q. Yang, Computational discretization algorithms for functional inequality
constrained optimization, Annals of Operations
Research, Vol. 98, 215/234 (2000)
8) M. Zibulevsky and M. Elad, L1-L2 optimizationin signal and image processing, IEEE Signal Processing Magazine, Vol. 27, 76/88 (2010)