目次
はじめに
スタンフォード大学には
機械学習を学ぶ上での第一歩として、
Introduction to Matrix Methods (EE103)という授業があります。
今回の記事では、
この授業の教科書である
Introduction to Applied Linear Algebraを
読んだ際の技術メモです。
この教科書は下記のリンクのページから
pdfをダウンロードすることができます。
本記事では、
上記の教科書の制約付き最小二乗法の部分のみのメモです。
他の部分に関しては、下記の記事を参照下さい。
制約付き最小二乗法
前回の記事では、通常の最小二乗法について説明しましたが、
今回は、線形の制約条件がついた場合の最小二乗法について説明します。
このような線形制約付きの最小二乗法も、
QR分解を使って、解くことができます。
制約付き最小二乗法は下記のように定式化されます。
Cは行列、dはベクトルです。
制約付き最小二乗法の解き方
この問題を解く、最もシンプルな方法は、
線形制約条件を目的関数とし、
下記の式のように、
複数目的関数における最小二乗法問題とし、
ラムダを非常に大きくした状態で最小二乗法を解くことで、
線形制約条件を満たした状態で、目的関数を最小化できます。
これをもう少し洗練された方法にしたのが、
ラグランジュの未定乗数法です。
先程の複数目的関数の最小二乗法問題の
目的関数をラグランジュ関数、
λをラグランジュ乗数と呼び、
元となった制約付き最小二乗問題の最適値x*と、
それに伴うラグランジュ乗数λ*は、
下記の条件を満たすことになります。
ここで、二つ目の項が0なのは、
ラグランジュ乗数でLを微分すると、
||Cx - d||^2 = 0となり、Cx = dで、
まさにはじめの線形制約を満たす条件となります。
１つ目の制約は、
二つ目の制約条件を満たした時に、
||Ax-b||^2の極値xだと言えるため、
その時のxは最適値になります。
ここで、ラグランジュ関数Lを展開して、
実際にxで微分すると、
となり、この式は下記のような線形方程式で表せます。
この式をKKT条件といい、
この線形方程式を解けば、最適値x*を求めることができます。
ちなみに左側の係数行列をKKT行列と呼びます。
線形方程式の解き方は下記を参照ください。
加えて、コスト関数が線形ではなく、
二次形式の場合も同じようにKKT条件を設定し、
KKT行列を用いた線形方程式を解けば、
解くことができます。
下記を参照ください。
制約付き最小二乗法の応用例
区分多項式によるフィッティング
状態空間を区分して、各区間に対してデータフィッティングをしたい場合、
制約付き最小二乗法が使えます。
例えば、下記のようにaという値を境に、
右と左で別々の多項式でデータフィッティングしたい場合、
aの地点において、yの値と方向(yの微分値)を同じにするために、
線形制約を置くことができます。
この線形制約付きの最小二乗法を解くことで、
上記のような区分された多項式フィッティングが可能になります。
このようなフィッティング手法を一般的にSplineといいます。
３次のスプライン補間に関しては、下記を参照ください。
では実際にJuliaを使って、
この区間多項式のフィッティングを実施してみたいと思います。
下記のサンプルコードでは、-1から1までのとある関数を
0で、左右２つに分けた区間多項式で近似します。
この時、左右の区間多項式が接する0でyの値と微分値が同じになるように
線形制約を加えます。
下記が先程の制約付き最小二乗法の各行列です。
Juliaを使って、この区間多項式フィッティングを解くコードは下記の通りです。
# # Piecewise Polynomial with constrainded least squares problem # # author: Atsushi Sakai # using PyPlot """ calc vandermonde matrix for polynomial calculation t: elements n: degree """ function calc_vandermonde_matrix(t,n) m = length(t) V = zeros(m,n) for i = 1:m, j = 1:n V[i,j] = t[i]^(j-1) end return V end """ solve constraint least squares min ||Ax - b||^2 s.t Cx = d """ function solve_constrained_least_squares(A,b,C,d) m, n = size(A) p, n = size(C) G = A'*A # Gram matrix KKT = [ 2*G C'; C zeros(p,p)] # KKT matrix xzhat = KKT \ [2*A'*b; d] return xzhat[1:n,:] end function main() println(PROGRAM_FILE," start!!") M = 70 # number of data point of one class N = 2*M # total number of data n = 4 # polynomial degree = n-1 Npl = 200 # number of polynomial samples xleft = rand(M) .- 1 xright = rand(M) x = [xleft; xright] y = x.^3 - x +0.4./(1.+25*x.^2) # set up matrix for constrained least square A = [ calc_vandermonde_matrix(xleft,n) zeros(M,n); zeros(M,n) calc_vandermonde_matrix(xright,n)] # display(A) b = y # println(b) C = [1 zeros(1,n-1) -1 zeros(1,n-1); 0 1 zeros(1,n-2) 0 -1 zeros(1,n-2)]; # println(C) d = zeros(2); # println(d) theta = solve_constrained_least_squares(A,b,C,d) xpl_left = linspace(-1, 0, Npl) ypl_left = calc_vandermonde_matrix(xpl_left, 4)*theta[1:n] xpl_right = linspace(0, 1, Npl) ypl_right = calc_vandermonde_matrix(xpl_right, 4)*theta[n+1:end] plot(x, y, ".k") plot(xpl_left, ypl_left, "-g") plot(xpl_right, ypl_right, "-r") axis("equal") show() println(PROGRAM_FILE," Done!!") end main()
上記を実行すると、
下記のような区間多項式フィッティングの近似結果が表示されます。
予算配分問題
予算の総数が決まっている状態で、
ある目標値にできるだけ近づける方法として、
制約付き最小二乗法を使うことができます。
投資における予算配分問題は、
ポートフォリオ最適化と呼ばれ、
制約付き最小二乗法の最も有名な応用例の一つになっています。
線形二次制御 (Linear Quadratic Control)
線形二次制御は制約付き最小二乗法を使った制御手法です。
以前、紹介したLQRは無限時間での最適制御でしたが、
今回の線形二次制御では有限時間での最適制御を実施します。
以前紹介したMPCの不等式制約が無い制御手法だとともいえます。
まず、下記のような線形運動モデルと観測モデルを考えます。
ここで状態量xを0にし、かつ、入力uをできるだけ小さくする制御を実施するために、
下記のような最適制御問題を考えます。
まずコスト関数としては、
J_outputという観測値(システムの出力)の二乗和と
J_input という入力値の二乗和を最小化する
複数の目的関数を設定します。
ρはこの２つの項のバランスを決定する重み値です。
制約条件に関しては、
先程の線形運動方程式と、
初期状態、終端状態の制約条件になります。
このように、制約条件は線形、
コスト関数は二次形式の形をしているため、
この最適制御手法は、線形二次制御と呼ばれています。
続いて、これらを制約付き最小二乗法として解きます。
まず、状態xと入力uを連結した拡張ベクトルを作成します。
そして、コスト関数に関しては、下記の行列を作成し、
制約条件に関しては、下記の２つの行列を作成すると、
先程の最適化問題は下記のように、
制約付き最小二乗法として定式化することができます。
yの値を0にするように制御したい場合は、
上記のbは0ベクトルになりますが、
ある目標値に制御したい場合は、bに目標値を設定することで、
任意の目標値に向けた制御が可能になります。
この制約付き最小二乗法は、
先程のKKT条件の線形方程式を解くことで解くことが可能です。
あとは、最適拡張ベクトルの値の入力値を使って、
システムの制御を実施することができます。
下記は実際にある線形モデルの制御において、
u = 0で制御したのも(Open loop制御)
と前述のLinear Quadratic Controlを使って制御を実施するJuliaコードです。
# # Linear quadratic control # using PyPlot function cls_solve(A,b,C,d) m, n = size(A) p, n = size(C) Q, R = qr([A; C]) Q = Matrix(Q) Q1 = Q[1:m,:] Q2 = Q[m+1:m+p,:] Qtil, Rtil = qr(Q2') Qtil = Matrix(Qtil) w = Rtil \ (2*Qtil'*Q1'*b - 2*(Rtil'\d)) return R \ (Q1'*b - Q2'*w/2) end function lqc(A, B, C, x_init, x_des, T, rho) n = size(A,1) m = size(B,2) p = size(C,1) q = size(x_init,2) Atil = [ kron(eye(T), C) zeros(p*T,m*(T-1)) ; zeros(m*(T-1), n*T) sqrt(rho)*eye(m*(T-1)) ] btil = zeros(p*T + m*(T-1), q) # We’ll construct Ctilde bit by bit Ctil11 = [ kron(eye(T-1), A) zeros(n*(T-1),n) ] - [ zeros(n*(T-1), n) eye(n*(T-1)) ] Ctil12 = kron(eye(T-1), B) Ctil21 = [eye(n) zeros(n,n*(T-1)); zeros(n,n*(T-1)) eye(n)] Ctil22 = zeros(2*n,m*(T-1)) Ctil = [Ctil11 Ctil12; Ctil21 Ctil22] dtil = [zeros(n*(T-1), q); x_init; x_des] z = cls_solve(Atil,btil,Ctil,dtil) x = [z[(i-1)*n+1:i*n,:] for i=1:T] u = [z[n*T+(i-1)*m+1 : n*T+i*m, :] for i=1:T-1] y = [C*xt for xt in x] return x, u, y end function main() println(PROGRAM_FILE," start!!") A = [ 0.855 1.161 0.667; 0.015 1.073 0.053; -0.084 0.059 1.022 ] B = [-0.076; -0.139; 0.342 ] C = [ 0.218 -3.597 -1.683 ] n = 3 p = 1 m = 1 x_init = [0.496; -0.745; 1.394] x_des = zeros(n,1) # open loop control u = 0 T = 100 yol = zeros(T,1) Xol = [ x_init zeros(n, T-1) ] for k=1:T-1 Xol[:,k+1] = A*Xol[:,k]; end yol = C*Xol # linear quadratic control rho = 0.2 x, u, y = lqc(A,B,C,x_init,x_des,T,rho) plot(yol', label="Open loop") plot(vcat(y...), label="Linear quadratic control") legend() show() println(PROGRAM_FILE," Done!!") end @time main()
上記のコードを実行すると、
下記のような制御結果が表示されます。
Open loopの制御は最後の制御値が0になっていませんが、
Linear Quadratic Controlを使って制御を実施したばあい、
最終値が0に制御されていることがわかります。
線形二次状態推定 (Linear Quadratic State Estimation)
前述の線形二次制御と同じ考え方で、
ある時系列の観測値から、制約付き最小二乗法を使って、
時系列の状態量を推定することができます。
これが線形二次状態推定です。
これはGPSの軌跡を計算する後処理などに利用されています。
同様の考え方を毎時刻毎に実行するのがカルマンフィルタです。
先程の同じ線形運動モデルと観測モデルを元に、
線形二次状態推定では下記のような最適化式を解きます。
ここで、コスト関数の１つ目の項は
観測値yと状態との差を表す、観測誤差vを最小化する項、
二つ目の項は、状態ベクトルの誤差を最小化する項です。
この最小化式を、先程のように解くことで、
時系列の状態量を推定することができます。
この線形二次状態推定においても、
データフィッティングの時と同様に、
最適用のデータとテスト用のデータを分けて、
クロスバリデーションをすることで、重みパラメータである
λを決定することが多いようです。
では実際にJuliaで線形二次状態推定を実施してみたいと思います。
シンプルなdouble integratorの運動モデルにおいて、
受け取った観測値から状態軌跡を推定するコードは下記のようになります。
# # Linear quadratic state estimation # # author: Atsushi Sakai # using PyPlot """ Get measurement """ function get_measurements(T) x = linspace(0,5,T) y = sin.(x).+rand(length(x))*0.5 m = [x y]' return m end """ solve constraint least squares min ||Ax - b||^2 s.t Cx = d """ function solve_constrained_least_squares(A,b,C,d) m, n = size(A) p, n = size(C) G = A'*A # Gram matrix KKT = [ 2*G C'; C zeros(p,p)] # KKT matrix xzhat = KKT \ [2*A'*b; d] return xzhat[1:n,:] end """ least square state estimation """ function least_square_state_estimation(A,B,C,y,T,lambda) n = size(A,1) m = size(B,2) p = size(C,1) Atil = [ kron(eye(T), C) zeros(T*p, m*(T-1)); zeros(m*(T-1), n*T) sqrt(lambda)*eye(m*(T-1)) ] btil = [ vcat(y...) ; zeros((T-1)*m) ] Ctil = [ ([ kron(eye(T-1), A) zeros(n*(T-1), n) ] + [ zeros(n*(T-1), n) -eye(n*(T-1)) ]) kron(eye(T-1), B) ] dtil = zeros(n*(T-1)) z = solve_constrained_least_squares(Atil, btil, Ctil, dtil) x = [ z[(i-1)*n+1:i*n] for i=1:T ] u = [ z[n*T+(i-1)*m+1 : n*T+i*m] for i=1:T-1 ] y = [ C*xt for xt in x ] # convert list to array x = hcat(x...) u = hcat(u...) y = hcat(y...) return x, u, y end function main() println(PROGRAM_FILE," start!!") T = 100 lambda = 10^3 A = [ eye(2) eye(2); zeros(2,2) eye(2) ] B = [ zeros(2,2); eye(2) ] C = [ eye(2) zeros(2,2) ] y = get_measurements(T) xest, uest, yest = least_square_state_estimation(A, B, C, y, T, lambda) plot(y[1,:],y[2,:],".k") plot(yest[1,:],yest[2,:],".-r") axis("equal") show() println(PROGRAM_FILE," Done!!") end main()
上記を実行すると、下記のような観測値と、
推定された軌跡が表示されます。
GItHubリポジトリ
前述のJuliaのコードは下記のリポジトリで公開しています。
参考資料
MyEnigma Supporters
もしこの記事が参考になり、
ブログをサポートしたいと思われた方は、
こちらからよろしくお願いします。
