MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

スタンフォード大学の学生が学ぶ制約付き最小二乗法とその応用

目次

 

はじめに

スタンフォード大学には

機械学習を学ぶ上での第一歩として、

Introduction to Matrix Methods (EE103)という授業があります。

 

今回の記事では、

この授業の教科書である

Introduction to Applied Linear Algebraを

読んだ際の技術メモです。

 

この教科書は下記のリンクのページから

pdfをダウンロードすることができます。

 

本記事では、

上記の教科書の制約付き最小二乗法の部分のみのメモです。

他の部分に関しては、下記の記事を参照下さい。

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

制約付き最小二乗法

前回の記事では、通常の最小二乗法について説明しましたが、

今回は、線形の制約条件がついた場合の最小二乗法について説明します。

このような線形制約付きの最小二乗法も、

QR分解を使って、解くことができます。

 

制約付き最小二乗法は下記のように定式化されます。

f:id:meison_amsl:20180813091200p:plain

Cは行列、dはベクトルです。

 

制約付き最小二乗法の解き方

この問題を解く、最もシンプルな方法は、

線形制約条件を目的関数とし、

下記の式のように、

複数目的関数における最小二乗法問題とし、

f:id:meison_amsl:20180813092258p:plain

ラムダを非常に大きくした状態で最小二乗法を解くことで、

線形制約条件を満たした状態で、目的関数を最小化できます。

 

これをもう少し洗練された方法にしたのが、

ラグランジュの未定乗数法です。

myenigma.hatenablog.com

 

先程の複数目的関数の最小二乗法問題の

目的関数をラグランジュ関数、

λをラグランジュ乗数と呼び、

元となった制約付き最小二乗問題の最適値x*と、

それに伴うラグランジュ乗数λ*は、

下記の条件を満たすことになります。

f:id:meison_amsl:20180816134603p:plain

 

ここで、二つ目の項が0なのは、

ラグランジュ乗数でLを微分すると、

||Cx - d||^2 = 0となり、Cx = dで、

まさにはじめの線形制約を満たす条件となります。

 

1つ目の制約は、

二つ目の制約条件を満たした時に、

||Ax-b||^2の極値xだと言えるため、

その時のxは最適値になります。

 

ここで、ラグランジュ関数Lを展開して、

実際にxで微分すると、

f:id:meison_amsl:20180816135915p:plain

となり、この式は下記のような線形方程式で表せます。

f:id:meison_amsl:20180816140555p:plain

この式をKKT条件といい、

この線形方程式を解けば、最適値x*を求めることができます。

ちなみに左側の係数行列をKKT行列と呼びます。

線形方程式の解き方は下記を参照ください。

myenigma.hatenablog.com

 

加えて、コスト関数が線形ではなく、

二次形式の場合も同じようにKKT条件を設定し、

KKT行列を用いた線形方程式を解けば、

解くことができます。

下記を参照ください。

myenigma.hatenablog.com

 

制約付き最小二乗法の応用例

区分多項式によるフィッティング

状態空間を区分して、各区間に対してデータフィッティングをしたい場合、

制約付き最小二乗法が使えます。

例えば、下記のようにaという値を境に、

右と左で別々の多項式でデータフィッティングしたい場合、

aの地点において、yの値と方向(yの微分値)を同じにするために、

線形制約を置くことができます。

f:id:meison_amsl:20180816125514p:plain

この線形制約付きの最小二乗法を解くことで、

上記のような区分された多項式フィッティングが可能になります。

このようなフィッティング手法を一般的にSplineといいます。

3次のスプライン補間に関しては、下記を参照ください。

myenigma.hatenablog.com

 

では実際にJuliaを使って、

この区間多項式のフィッティングを実施してみたいと思います。

下記のサンプルコードでは、-1から1までのとある関数を

0で、左右2つに分けた区間多項式で近似します。

この時、左右の区間多項式が接する0でyの値と微分値が同じになるように

線形制約を加えます。

下記が先程の制約付き最小二乗法の各行列です。

f:id:meison_amsl:20180825211031p:plain

f:id:meison_amsl:20180825211041p:plain

  

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()

上記を実行すると、

下記のような区間多項式フィッティングの近似結果が表示されます。

f:id:meison_amsl:20180825210831p:plain

 

予算配分問題

予算の総数が決まっている状態で、

ある目標値にできるだけ近づける方法として、

制約付き最小二乗法を使うことができます。

 

投資における予算配分問題は、

ポートフォリオ最適化と呼ばれ、

制約付き最小二乗法の最も有名な応用例の一つになっています。

myenigma.hatenablog.com

 

線形二次制御 (Linear Quadratic Control)

線形二次制御は制約付き最小二乗法を使った制御手法です。

以前、紹介したLQRは無限時間での最適制御でしたが、

myenigma.hatenablog.com

今回の線形二次制御では有限時間での最適制御を実施します。

以前紹介したMPCの不等式制約が無い制御手法だとともいえます。

myenigma.hatenablog.com

 

まず、下記のような線形運動モデルと観測モデルを考えます。

f:id:meison_amsl:20180825103224p:plain

ここで状態量xを0にし、かつ、入力uをできるだけ小さくする制御を実施するために、

下記のような最適制御問題を考えます。

f:id:meison_amsl:20180825103618p:plain

まずコスト関数としては、

J_outputという観測値(システムの出力)の二乗和と

J_input という入力値の二乗和を最小化する

複数の目的関数を設定します。

ρはこの2つの項のバランスを決定する重み値です。

制約条件に関しては、

先程の線形運動方程式と、

初期状態、終端状態の制約条件になります。

このように、制約条件は線形、

コスト関数は二次形式の形をしているため、

この最適制御手法は、線形二次制御と呼ばれています。

 

続いて、これらを制約付き最小二乗法として解きます。

まず、状態xと入力uを連結した拡張ベクトルを作成します。

f:id:meison_amsl:20180825105832p:plain

そして、コスト関数に関しては、下記の行列を作成し、

f:id:meison_amsl:20180825105151p:plain

制約条件に関しては、下記の2つの行列を作成すると、

f:id:meison_amsl:20180825105232p:plain

先程の最適化問題は下記のように、

制約付き最小二乗法として定式化することができます。

f:id:meison_amsl:20180825105358p:plain

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()

上記のコードを実行すると、

下記のような制御結果が表示されます。

f:id:meison_amsl:20180901095155p:plain

Open loopの制御は最後の制御値が0になっていませんが、

Linear Quadratic Controlを使って制御を実施したばあい、

最終値が0に制御されていることがわかります。

 

線形二次状態推定 (Linear Quadratic State Estimation)

前述の線形二次制御と同じ考え方で、

ある時系列の観測値から、制約付き最小二乗法を使って、

時系列の状態量を推定することができます。

これが線形二次状態推定です。

これはGPSの軌跡を計算する後処理などに利用されています。

同様の考え方を毎時刻毎に実行するのがカルマンフィルタです。

myenigma.hatenablog.com

 

先程の同じ線形運動モデルと観測モデルを元に、

線形二次状態推定では下記のような最適化式を解きます。

f:id:meison_amsl:20180825112017p:plain

ここで、コスト関数の1つ目の項は

観測値yと状態との差を表す、観測誤差vを最小化する項、

f:id:meison_amsl:20180825112131p:plain

二つ目の項は、状態ベクトルの誤差を最小化する項です。

f:id:meison_amsl:20180825112403p:plain

この最小化式を、先程のように解くことで、

時系列の状態量を推定することができます。

 

この線形二次状態推定においても、

データフィッティングの時と同様に、

最適用のデータとテスト用のデータを分けて、

クロスバリデーションをすることで、重みパラメータである

λを決定することが多いようです。

 

では実際に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()

上記を実行すると、下記のような観測値と、

推定された軌跡が表示されます。

f:id:meison_amsl:20180826095721p:plain

 

GItHubリポジトリ

前述のJuliaのコードは下記のリポジトリで公開しています。

github.com

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

 

MyEnigma Supporters

もしこの記事が参考になり、

ブログをサポートしたいと思われた方は、

こちらからよろしくお願いします。

myenigma.hatenablog.com