notebooks/RL_LQR_estimation_and_robustness.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LQRにおけるサンプル効率\n",
    "\n",
    "参考：\n",
    "* [On the Sample Complexity of the Linear Quadratic Regulator](https://arxiv.org/abs/1710.01688)\n",
    "\n",
    "モデルのダイナミクス行列$A$と$B$がわかっていればLQRでダイナミクスを解くことができますが，一般にはそうではありません．\n",
    "わからないときは，ダイナミクスを推定したあとに最適制御を実行する必要があります．\n",
    "\n",
    "表記：\n",
    "$$\n",
    "\\begin{array}{ll}\n",
    "\\operatorname{minimize} & \\mathbb{E}\\left[\\frac{1}{T} \\sum_{t=1}^T x_t^* Q x_t+u_{t-1}^* R u_{t-1}\\right] \\\\\n",
    "\\text { subject to } & x_{t+1}=A x_t+B u_t+w_t\n",
    "\\end{array}\n",
    "$$\n",
    "* $Q$は$n\\times n$の行列\n",
    "* $w_t\\in \\mathbb{R}^n$は平均０で分散$\\Sigma_w$のガウシアンノイズ行列\n",
    "* $G_T=\\left[\\begin{array}{llll}A^{T-1} B & A^{T-2} B & \\ldots & B\\end{array}\\right] \\quad$ and $\\quad F_T=\\left[\\begin{array}{llll}A^{T-1} & A^{T-2} & \\ldots & I_n\\end{array}\\right]$とする．\n",
    "\n",
    "\n",
    "## ダイナミクスの推定\n",
    "\n",
    "ガウシアンノイズによってシステムをロールアウトして，行列$A$と$B$を推定することを考えます．つまり，\n",
    "\n",
    "$$\n",
    "x_{t+1}^{(\\ell)}=A x_t^{(\\ell)}+B u_t^{(\\ell)}+w_t^{(\\ell)} \\text { with } w_t^{(\\ell)} \\stackrel{\\text { i.i.d. }}{\\sim} \\mathcal{N}\\left(0, \\sigma_w^2 I_n\\right) \\text { and } u_t^{(\\ell)} \\stackrel{\\text { i.i.j. }}{\\sim} \\mathcal{N}\\left(0, \\sigma_u^2 I_p\\right) \\text {. }\n",
    "$$\n",
    "\n",
    "として，データセット$\\left\\{\\left(x_t^{(\\ell)}, u_t^{(\\ell)}\\right): 1 \\leq \\ell \\leq N, 0 \\leq t \\leq T\\right\\}$を集め，システムを\n",
    "$$\n",
    "(\\widehat{A}, \\widehat{B}) \\in \\arg \\min _{(A, B)} \\sum_{\\ell=1}^N \\frac{1}{2}\\left\\|A x_{T-1}^{(\\ell)}+B u_{T-1}^{(\\ell)}-x_T^{(\\ell)}\\right\\|_2^2\n",
    "$$\n",
    "として推定します．（これは解析を簡単にするために，ロールアウトで一個のサンプルが得られる状況だけ考えます．一方で，ロールアウトの全てのサンプルを使用する場合はBootstrap法として論文中に記載があるので，そちらを参照してください）\n",
    "\n",
    "このとき，$N \\geq 8(n+p)+16 \\log (4 / \\delta)$ならば確率$1-\\delta$以上で\n",
    "$$\n",
    "\\|\\widehat{A}-A\\|_2 \\leq \\frac{16 \\sigma_w}{\\sqrt{\\lambda_{\\min }\\left(\\sigma_u^2 G_T G_T^*+\\sigma_w^2 F_T F_T^*\\right)}} \\sqrt{\\frac{(n+2 p) \\log (36 / \\delta)}{N}}\n",
    "$$\n",
    "であり，また，\n",
    "$$\n",
    "\\|\\widehat{B}-B\\|_2 \\leq \\frac{16 \\sigma_w}{\\sigma_u} \\sqrt{\\frac{(n+2 p) \\log (36 / \\delta)}{N}}\n",
    "$$\n",
    "が成り立つことを示しましょう．\n",
    "\n",
    "---\n",
    "\n",
    "**別の式表記**\n",
    "\n",
    "証明の前に，式展開を簡単にするために表記を導入します．\n",
    "\n",
    "$\n",
    "\\Theta:=\\left[\\begin{array}{ll}\n",
    "A & B\n",
    "\\end{array}\\right]^* \\in \\mathbb{R}^{(n+p) \\times n}\n",
    "$\n",
    "および\n",
    "$z_t:=\\left[\\begin{array}{l}x_t \\\\ u_t\\end{array}\\right] \\in \\mathbb{R}^{n+p}$\n",
    "とします．\n",
    "\n",
    "このとき，\n",
    "\n",
    "$$\n",
    "x_{t+1}^*=z_t^* \\Theta+w_t^*\n",
    "$$\n",
    "\n",
    "が成立し，一度のロールアウトでは\n",
    "\n",
    "$$\n",
    "X:=\\left[\\begin{array}{c}\n",
    "x_1^* \\\\\n",
    "x_2^* \\\\\n",
    "\\vdots \\\\\n",
    "x_T^*\n",
    "\\end{array}\\right], Z:=\\left[\\begin{array}{c}\n",
    "z_0^* \\\\\n",
    "z_1^* \\\\\n",
    "\\vdots \\\\\n",
    "z_{T-1}^*\n",
    "\\end{array}\\right], W:=\\left[\\begin{array}{c}\n",
    "w_0^* \\\\\n",
    "w_1^* \\\\\n",
    "\\vdots \\\\\n",
    "w_{T-1}^*\n",
    "\\end{array}\\right]\n",
    "$$\n",
    "\n",
    "が得られます．\n",
    "\n",
    "このときのシステムのダイナミクスは$X=Z\\Theta + W$で与えられることに注意して，ロールアウトのデータセットを$X_N, Z_N, W_N$とすれば，（ここで$X_N=\\left[\\begin{array}{llll}x_T^{(1)} & x_T^{(2)} & \\ldots & x_T^{(N)}\\end{array}\\right]^*$）\n",
    "\n",
    "$$X_N=Z_N \\Theta+W_N$$\n",
    "\n",
    "で与えれれます．よって，$Z_N^*Z_N$がinvertibleであるとすれば，パラメータの推定は\n",
    "\n",
    "$$\n",
    "\\widehat{\\Theta}=\\left(Z_N^* Z_N\\right)^{-1} Z_N^* X_N=\\Theta+\\left(Z_N^* Z_N\\right)^{-1} Z_N^* W_N\n",
    "$$\n",
    "\n",
    "であり，推定誤差は\n",
    "\n",
    "$$\n",
    "E:=\\widehat{\\Theta}-\\Theta=\\left(Z_N^* Z_N\\right)^{-1} Z_N^* W_N\n",
    "$$\n",
    "\n",
    "になります．\n",
    "\n",
    "---\n",
    "\n",
    "$G_T$と$F_T$の定義を思い出しましょう．$G_T=\\left[\\begin{array}{llll}A^{T-1} B & A^{T-2} B & \\ldots & B\\end{array}\\right], \\quad F_T=\\left[\\begin{array}{llll}A^{T-1} & A^{T-2} & \\ldots & I_n\\end{array}\\right]$だったので，\n",
    "\n",
    "$$\n",
    "x_T=G_T\\left[\\begin{array}{c}\n",
    "u_0 \\\\\n",
    "u_1 \\\\\n",
    "\\vdots \\\\\n",
    "u_{T-1}\n",
    "\\end{array}\\right]+F_T\\left[\\begin{array}{c}\n",
    "w_0 \\\\\n",
    "w_1 \\\\\n",
    "\\vdots \\\\\n",
    "w_{T-1}\n",
    "\\end{array}\\right]\n",
    "$$\n",
    "\n",
    "が成り立ちます．$u_t \\sim \\mathcal{N}\\left(0, \\sigma_u^2 I_p\\right)$だったので，\n",
    "\n",
    "$$\n",
    "\\left[\\begin{array}{l}\n",
    "x_T \\\\\n",
    "u_T\n",
    "\\end{array}\\right] \\sim \\mathcal{N}\\left(0,\\left[\\begin{array}{cc}\n",
    "\\sigma_u^2 G_T G_T^*+\\sigma_w^2 F_T F_T^* & 0 \\\\\n",
    "0 & \\sigma_u^2 I_p\n",
    "\\end{array}\\right]\\right)\n",
    "$$\n",
    "\n",
    "です．$F_T F_T^* \\succ 0$なので，$\\sigma_u$と$\\sigma_w$が正であれば，これは退化していない分布になります．\n",
    "\n",
    "これを踏まえて，次の補題を使います．\n",
    "\n",
    "---\n",
    "\n",
    "$\\delta \\in (0, 1)$かつ$N \\geq 2 \\log(1/\\delta)$とする．\n",
    "$f_k \\in \\mathbb{R}^m, g_k \\in \\mathbb{R}^n$を$f_k \\sim \\mathcal{N}\\left(0, \\Sigma_f\\right)$ and $g_k \\sim \\mathcal{N}\\left(0, \\Sigma_g\\right)$なる独立確率ベクトルとする．\n",
    "このとき，確率$1-\\delta$以上で，\n",
    "\n",
    "$$\n",
    "\\left\\|\\sum_{k=1}^N f_k g_k^*\\right\\|_2 \\leq 4\\left\\|\\Sigma_f\\right\\|_2^{1 / 2}\\left\\|\\Sigma_g\\right\\|_2^{1 / 2} \\sqrt{N(m+n) \\log (9 / \\delta)}\n",
    "$$\n",
    "\n",
    "が成立する．\n",
    "\n",
    "証明はBernsteinをいじるとでてくるよ．Appendix A参照．\n",
    "\n",
    "---\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "この補題を使うと，$X$が$n_1\\times N$の$\\mathcal{N}(0, 1)$に従う確率分布， $Y$が$N\\times n_2$の$\\mathcal{N}(0, 1)$に従う確率分布とすると，確率$1-\\delta$以上で，\n",
    "\n",
    "$$\n",
    "\\|X Y\\|_2 \\leq 4 \\sqrt{N\\left(n_1+n_2\\right) \\log (9 / \\delta)}\n",
    "$$\n",
    "\n",
    "が成り立つ．続いて，次の行列についての特異値のバウンドも便利です．\n",
    "\n",
    "---\n",
    "\n",
    "$X \\in \\mathbb{R}^{N \\times n}$を$\\mathcal{N}(0, 1)$に従う独立確率行列とすれば，確率$1-\\delta$以上で，\n",
    "\n",
    "$$\n",
    "\\sqrt{\\lambda_{\\min }\\left(X^* X\\right)} \\geq \\sqrt{N}-\\sqrt{n}-\\sqrt{2 \\log (1 / \\delta)}\n",
    "$$\n",
    "\n",
    "が成り立つ．\n",
    "\n",
    "---\n",
    "\n",
    "この２つを合体すると，次がわかります．\n",
    "\n",
    "---\n",
    "\n",
    "* $z_1, \\dots, z_N \\in \\mathbb{R}^n$を$\\mathcal{N}(0, \\Sigma)$に従うとして，$\\Sigma$をinvertibleとする．\n",
    "* $Z^*:=\\left[\\begin{array}{lll}z_1 & \\ldots & z_N\\end{array}\\right]$とする．\n",
    "* $W\\in \\mathbb{R}^{N\\times p}$の$\\mathcal{N}(0, \\sigma_w^2)$に従うとして，$Z$と独立とする．\n",
    "* $E:=\\left(Z^* Z\\right)^{\\dagger} Z^* W$とする．\n",
    "\n",
    "\n",
    "このとき，$N \\geq 8 n+16 \\log (2 / \\delta)$のとき，任意の行列$Q$に対して，確率$1-\\delta$以上で\n",
    "\n",
    "$$\n",
    "\\|Q E\\|_2 \\leq 16 \\sigma_w\\left\\|Q \\Sigma^{-1 / 2}\\right\\|_2 \\sqrt{\\frac{(n+p) \\log (18 / \\delta)}{N}}\n",
    "$$\n",
    "\n",
    "が成り立つ．\n",
    "\n",
    "証明は論文参照．\n",
    "\n",
    "---\n",
    "\n",
    "以上を使えば，$\\hat{A}-A$および$\\hat{B}-B$の推定誤差のバウンドは簡単に出せます．\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "---\n",
    "\n",
    "**バウンドのコメント**\n",
    "\n",
    "$$\n",
    "\\|\\widehat{A}-A\\|_2 \\leq \\frac{16 \\sigma_w}{\\sqrt{\\lambda_{\\min }\\left(\\sigma_u^2 G_T G_T^*+\\sigma_w^2 F_T F_T^*\\right)}} \\sqrt{\\frac{(n+2 p) \\log (36 / \\delta)}{N}}\n",
    "$$\n",
    "$$\n",
    "\\|\\widehat{B}-B\\|_2 \\leq \\frac{16 \\sigma_w}{\\sigma_u} \\sqrt{\\frac{(n+2 p) \\log (36 / \\delta)}{N}}\n",
    "$$\n",
    "\n",
    "* 学習すべきパラメータは$n(n+p)$個ありますが，今回のバウンドは$(n+2p)$の依存なので，$n$と$p$についてtightです．\n",
    "* $\\lambda_{\\min }\\left(\\sigma_u^2 G_T G_T^*+\\sigma_w^2 F_T F_T^*\\right)$は入力やノイズが状態全体にどれくらい増大されるか？の量になります．"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 推定したダイナミクスでの制御\n",
    "\n",
    "ロバストな制御を実現するためにSystem Level Synthesisを使います．つまり，状態と入力をノイズだけで次のように表します．\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "& x_k=\\sum_{t=1}^k(A+B K)^{k-t} w_{t-1}, \\\\\n",
    "& u_k=\\sum_{t=1}^k K(A+B K)^{k-t} w_{t-1}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "すると，$\\Phi_x(k):=(A+B K)^{k-1}$および$\\Phi_u(k):=K(A+B K)^{k-1}$とすれば，\n",
    "\n",
    "$$\\left[\\begin{array}{l}x_k \\\\ u_k\\end{array}\\right]=\\sum_{t=1}^k\\left[\\begin{array}{l}\\Phi_x(k-t+1) \\\\ \\Phi_u(k-t+1)\\end{array}\\right] w_{t-1}$$\n",
    "\n",
    "が成り立ちます．$\\left\\{\\Phi_x(k), \\Phi_u(k)\\right\\}$は$K$による閉ループシステムのシステム応答と呼ばれます．\n",
    "この$\\Phi$について\n",
    "\n",
    "$$\n",
    "\\Phi_x(k+1)=A \\Phi_x(k)+B \\Phi_u(k), \\Phi_x(1)=I, \\quad \\forall k \\geq 1\n",
    "$$\n",
    "\n",
    "が成り立ちます．\n",
    "\n",
    "また，無限ホライゾンの場合はz変換を使って表記すると便利です．$\\boldsymbol{\\Phi}_x(z)=\\sum_{k=1}^{\\infty} \\Phi_x(k) z^{-k}$とすると，上の制約は\n",
    "\n",
    "$$\n",
    "\\left[\\begin{array}{ll}\n",
    "z I-A & -B\n",
    "\\end{array}\\right]\\left[\\begin{array}{l}\n",
    "\\mathbf{\\Phi}_x \\\\\n",
    "\\mathbf{\\Phi}_u\n",
    "\\end{array}\\right]=I\n",
    "$$\n",
    "\n",
    "として表現できます．このときのゲインは$\\mathbf{u}=\\mathbf{K} \\mathbf{x}$ is given by $\\mathbf{K}=\\boldsymbol{\\Phi}_u \\boldsymbol{\\Phi}_x^{-1}$です．\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}