HyperX Alloy FPS Pro

スイスでは3月中旬にCOVID19の感染爆発に伴うロックダウンが実施されました。私もWFHに入ったのですが、ノートパソコン・スタンドを使った環境となったので 外付けキーボードが必要となり、KingstonのHyperX FPS Proを購入しました。 使いはじめて4ヶ月ほど経ったのでレビューします。

購入の動機

  • スイスで手に入るキーボードのうち、QWERTZなCHキーボードではなく、
    QWERTYなキーボードがほしかった
  • メカニカルキーボードに興味あった
  • キーボードは卓上を占拠しないコンパクトなものが好き
  • 予算100CHF(約1万円以内) / 納期1週間以内

使ってみて

  • 購入の動機は達成。テンキーなしにも慣れたので、さらに小さいキーボードに移行できるかも。
  • 一番スタンダードな赤軸(Cherry)ということですが、やはりペコペコキーボードより音が大きい。慣れてきた結果、底打ち(※)はしなくなりましたが、 静かなオフィスで使うのは同僚に遠慮して躊躇するかも。
    ※ どうやら「キーを必要以上に押し込んでしまい、大きな接触音がすること」を言うらしい
  • 購入した製品に不具合があるのか、たまにキーが押下固定されてしまい、"dddddddd...."などと虚無を量産することがある。redditでも何件か同じような不具合のコンプレインを目にした。
    • 別のキーを押すと開放されるが、その後、そのキーはキーボードにを再接続するまで使えない
    • 再現条件不明。同じキーを連打している時に起きやすい。2週間に1回程度の頻度。
    • 最新ファームウェア(2.1.1.2) はこちらから. アップデートのご利益なし。
      ->HyperX Alloy FPS Pro Firmware Update / Alloy FPS Pro | HyperX
    • しばらくは我慢してましたが、キーの誤動作で我慢ならん自体が発生したので、サポートに連絡したら代替品を出荷してくれた。到着待ちです。
      誤動作の原因が電気信号っぽい気がする。代替品作戦、地球環境に悪いし付録ケーブル含めてノイズとか確認したほうが良いと思うんですがどうですか>kingstonさん
  • 付属のUSBケーブルが硬くて長め。自分でUSB-mini B用意したほうが良いかも。

逐次ベイズ推定とカルマンフィルタ

カルマンフィルタは逐次ベイズフィルタの特殊系として言われますが、式を眺めていても結びつきがさっぱりイメージできなかったので数式いじりをしました。 もっとしっかりした入門は他でたくさん見かけますので、 参考文献として記事のリンクしておきます。

数式操作に関しては須山さんの本を参考にしました。

集合・位相・圏 数学の言葉への最短コース
ベイズ推論による機械学習入門
  • 作者: 須山敦志
  • 出版社/メーカー: 講談社
  • 発売日: 2017/10/21

逐次的ベイズ推定から見たカルマン・フィルタ

セオリーと違ってカルマンフィルタで得たい出力値 状態量 $ \hat{x}_{k|k} $ の位置付けから入ります。理解を簡単にするため、システムとしてD次元の酔歩 $ x_k \in \mathbb{R}^D $と 酔歩する粒子の観測誤差を含む観測値 $ y_{k} \in \mathbb{R}^D $ で考察します。

\[ \begin{array}{rll} x_{n+1} &= x_{n} + w_n & \qquad w_n \sim \mathcal{N}(0, Q) \quad Q \in \mathbb{R}^{D \times D} \\ y_{n} &= x_{n} + v_n & \qquad v_n \sim \mathcal{N}(0, R) \quad R \in \mathbb{R}^{D \times D} \\ \end{array} \]
観測値 $ y_n $ は真の平均値のわからない正規分布 $ \mathcal{N}(x_n, R) $ から取り出されたもので、 このn=1のデータからは、せいぜい $ \hat{x}_n = y_n $ とするのが関の山です。 しかし、うまく $ x_{n+1} - x_{n} \sim \mathcal{N}(0, Q) $ と、 観測したデータ $ Y_n = \{Y_{n-1}, y_{n}\} $ を使って、より高い精度で $ x_n $ を推定することを考えます。

カルマンフィルタの式

システムをざっくり簡略化したので、カルマンフィルタの式もザクッと簡単になります。この式は逐次ベイズ推定から求めた結果とあとから比較するために掲載しました。

予測のステップ
\[ \begin{array}{rl} \hat{x}_{k|k-1} &= \hat{x}_{k-1|k-1} \\ P_{k|k-1} &= P_{k-1|k-1} + Q \\ \end{array} \]
フィルタリングのステップ
\[ \begin{array}{rl} e_k &= y_k - \hat{x}_{k|k-1} \\ S_k &= R + P_{k|k-1} \\ K_k &= P_k S_k^{-1} = P_k (R + P_{k|k-1} )^{-1}\\ \hat{x}_{k|k} &= \hat{x}_{k|k-1} + K_k e_k\\ P_{k|k} &= (I - K_k) P_{k|k-1} \\ \end{array} \]

フィルタステップ

観測データ $ Y_n $から状態量 $ x_n $ の確率密度分布 $ p(x_n|Y_n) $ を求めます。確率密度分布にはベイズの定理から下記の性質があります。赤字では $ y_n $ と $ Y_{n-1} $ の独立性(マルコフ性)を用いました。

\[ \begin{array}{rll} p(x_n \vert Y_n) & = p(x_n | y_n, Y_{n-1}) \\ & \propto \color{red}{p(y_n | x_n, Y_{n-1})} p(x_n | Y_{n-1})\\ & = \color{red}{p(y_n | x_n)} p(x_n | Y_{n-1})\\ \end{array} \]

ここで、事前分布として $ p(x_n | Y_{n-1}) = \mathcal{N}(x_n | \hat{x}_{n \vert n-1}, P_{n \vert n-1}) $ をカルマンフィルタでは仮定する。このとき、$ p(y_n \vert x_n) = \mathcal{N}(y_n | x_n, R) $ であるので、

\[ \begin{array}{rll} \ln \left(p(y_n | x_n) p(x_n | Y_{n-1}) \right) &= \ln \mathcal{N}(y_n | x_n, k) + \ln \mathcal{N}(x_n | \hat{x}_{n|n-1}, P_{n|n-1}) \\ &= - \frac{1}{2} \left\{ x_n^T (\color{green}{R^{-1}+P_{n|n-1}^{-1}}) x_n - 2 x_n^T \color{blue}{(R^{-1} y_n + P_{n|n-1}^{-1} \hat{x}_{n|n-1}}) \right\} + \text{const.} \end{array} \]
この数式の形状により $ p(x_n \vert Y_n) $ がガウス分布であることがわかるので (須山本の手法です)、
\[ p(x_n \vert Y_n) = \mathcal{N}(x_n \vert \hat{x}_{n|n}, P_{n|n}) \\ \]
とおくと
\[ \begin{array}{rl} \ln p(x_n \vert Y_n) &= \left( -\frac{1}{2}(x_n - \hat{x}_{n|n})^T P_{n|n}^{-1} (x_n - \hat{x}_{n|n}) - \ln | P_{n|n}| \right) + \text{const}\\ & = -\frac{1}{2} \left( x_n^T \color{green}{P_{n|n}^{-1}} x_n - 2 x_n^T \color{blue}{P_{n|n}^{-1} \hat{x}_{n|n}} \right) + \text{const.} \end{array} \]
ここから各パラメータの比較をします。すると
\[ \begin{array}{crl} &P^{-1}_{n|n} &= \color{green}{R^{-1} + P_{n|n-1}^{-1}}\\ \Leftrightarrow & P_{n|n} &= \left( 1- P_{n|n-1}(P_{n|n-1}+R)^{-1} \right) P_{n|n-1} \\ && = \left( 1- K_n \right) P_{n|n-1} \\ \end{array} \]
と誤差共分散の更新式が求まります。さらに目当ての推定値が
\[ \begin{array}{crl} & P_{n|n}^{-1} \hat{x}_{n|n} &= \color{blue}{ R^{-1} y_n + P_{n|n-1}^{-1} \hat{x}_{n|n-1} }\\ \Leftrightarrow & \hat{x}_{n|n} &= P_{n|n} \left( R^{-1} y_n + P_{n|n-1}^{-1} \hat{x}_{n|n-1} \right) \\ && = \hat{x}_{n|n-1} + P_{n|n-1}(P_{n|n-1}+R)^{-1} (y_n - \hat{x}_{n|n-1}) \\ && = \hat{x}_{n|n-1} + K_n e_n \\ \text{where} & K_n &= P_{n|n-1}(P_{n|n-1}+R)^{-1} \\ & e_n &= y_n - \hat{x}_{n|n-1} \end{array} \]

と求まり、$ K $が、このシステムにおける最適カルマンゲインになっているなど、カルマンフィルタの更新ステップの式になっています。 $ \hat{x}_{n|n} $は、条件付き確率密度分布 $ p(x_n \vert Y_n) $ に従う、$ x_n $ の条件付き期待値 $ \mathbb{E}(x_n|Y_n) $ ですので、 これ自体が最小分散推定になっています(証明は省略. 参考文献[4]参照)。

予測のステップ

先ほどのフィルタリングのステップでは、事前分布として$ p(x_n | Y_{n-1}) = \mathcal{N}(x_n | \hat{x}_{n \vert n-1}, P_{n \vert n-1}) $を用いていました。この事前分布のパラメータ、$ \hat{x}_{n|n-1} $と$ P_{n|n-1} $を求めます。 $ x_{n+1} - x_{n} \sim \mathcal{N}(0, Q) $ であることから、 その1ステップ前のフィルタリングステップで求めた粒子位置に関する確率密度分布 $ p(x_{n-1}|Y_{n-1}) = \mathcal{N}(x_{n-1}|\hat{x}_{n-1|n-1}, P_{n-1|n-1}) $ に対して、 正規分布の加法性から

\[ \begin{array}{rl} \hat{x}_{n|n-1} &= \hat{x}_{n-1|n-1}\\ P_{n|n-1} &= P_{n-1|n-1} + Q\\ \end{array} \]
となります。そのままですね。

初期値

さて、同じようにしてフィルタリングステップ→予測ステップ→フィルタリングと、ステップをどんどんさかのぼっていくと、最初の位置 $ \hat{x}_{0|0} $, そしてその位置に関する共分散 $ P_{0|0} $ に行き着きます。 これはハイパーパラメータというか、えいやと決めるものになります。

システムがもう少し複雑な場合

前節では酔歩を仮定しましたが、一般の線形システムとして、

\[ \begin{array}{rl} x_{n+1} &= A x_n + B u_n + w_n \\ y_n &= C x_n + v_n \end{array} \]
について考えます。このとき、
\[ \begin{array}{rl} p(x_{n+1}|x_{n}) &= \mathcal{N}(A x_n + B u_n, Q) \\ p(y_n|x_n) &= \mathcal{N}(C x_n, R) \end{array} \]
であって、他の確率分布についても正規分布の再生性をふまえて
\[ \begin{array}{rl} x &\sim \mathcal{N}(\mu, \Sigma) \\ \Rightarrow G x &\sim \mathcal{N}(G \mu, G^T \Sigma G) \end{array} \]
に注意すれば、もとのカルマンフィルタの式になるのは自明ですね。

参考文献

  1. シンプルなモデルとイラストでカルマンフィルタを直観的に理解してみる - Qiita
  2. ベイズの定理からカルマンフィルタを導出する - ssk tech blog
  3. データ同化|カルマンフィルタと尤度 - ari23の研究ノート
  4. 知能制御システム学 画像追跡 (3) ― ベイズ推定とパーティクルフィルタ ― 東北大学 大学院情報科学研究科 鏡 慎吾 准教授
  5. カルマンフィルタの誤差共分散の意味 - 夜間飛行

Numpy使いのためのEigenチートシート

Eigen

v3.4からスライスが便利になるそうだが、最新Stable版の利用を前提にする。

時期が来たらこっそりエントリー書き換えるかも。

git clone https://gitlab.com/libeigen/eigen.git
git checkout tags/3.3.7

導入

#include <Eigen/Core>
using namespace Eigen;

基本

ベクトル - 1D Array - Vector

  • numpyのデフォルトはfloat64(=double)らしいので、それに倣う。
  Numpy Eigen
配列確保 a = np.empty(size) VectorXd a(size);
ゼロベクトル a = np.zeros(size) VectorXd a = VectorXd::Zero(size);
1ベクトル a = np.ones(size) VectorXd a = VectorXd::Ones(size);
fill a.fill(1) a.fill(1);
ベクトル長 a.shape[0] a.size();
要素アクセス a[0]
a[-1]
a(0)
a(a.size()-1)
スライス・アクセス a[a:a+n]
a[:n]
a[-n:]
a.segment(a,n)
a.head(n)
a.tail(n)
要素べき乗 a**n
a**2
a.array().pow(n)
a.array().square()
要素ルート a**0.5 a.array().sqrt()
要素総和 np.sum(a) a.sum()

行列 - 2D Array - Matrix

  Numpy Eigen
配列確保 a = np.empty((r,c)) MatrixXd a(r,c);
ゼロ行列 a = np.zeros((r,c)) MatrixXd a = MatrixXd::Zero(r,c);
1s 行列 a = np.ones((r,c)) MatrixXd a = MatrixXd::Ones(r,c);
サイズ a.shape a.rows(), a.cols()
要素アクセス a[r,c] a(r,c)
スライス・アクセス a[r:r+n,c:c+m] a.block<n,m>(r,c)
a.block(r,c,n,m)
列アクセス a[r] a.row(r)
行アクセス a[:,c] a.col(c)
flatten#1 a.flatten('F') Map<const VectorXd>(a.data(), a.size())
flatten#2 a.flatten() MatrixXd b = a.transpose();
Map<const VectorXd>(b.transpose().data(), b.size());
reshape
- MatrixXdのみ
a.reshape((r,c)) a.resize(r,c)
// 非破壊でreshapeするなら a.reshaped(r,c)
reshape a.reshape((r,c)) MatrixXd x = Map<Matrix<double, r, c>>(a.data());
transepose a.T a.transpose()
逆行列 np.linalg.inv(a) a.inverse()

配列構築・スニペット系

Eigenにはないnumpy/scipyの便利関数たち

  Numpy Eigen
vstack np.vstack([a, b])
MatrixXd c(a.rows()+b.rows(), a.cols());
c << a, b;
hstack np.hstack([a, b])
MatrixXd c(a.rows(), a.cols()+b.cols());
c << a, b;
対角行列 np.diag([0, 1, 2, 3])
MatrixXd v(4,1);
v << 0.,1.,2.,3;
MatrixXd m = v.asDiagonal();
block Diag # Scipy # numpycpp参照

参考

SympyでC++/Eigenな運動方程式の生成

ちょっと複雑な運動方程式をCやPythonコード(numba)に落とすとき、Sympyが便利です。Sympyのコード生成にちょっと癖があったのでメモ代わりに残します。 例として倒立振子をモデルにした。

モデル

欲しい式

非線形な運動方程式に対して、

\[ \dot{q} = f(q, u) \]

テーラー展開で一次線形近似した式の赤いところと青いところを計算してくれる式が欲しい。

\[ \dot{q} = \color{red}{f(q_0, u_0)} + \color{blue}{\frac{\partial f}{\partial q} \Big\vert_{q_0 u_0}} q +\color{blue}{ \frac{\partial f}{\partial u}\Big\vert_{q_0 u_0}} u \]

以降簡単のため、

\[ \begin{array}{rl} A_c(q_0, u_0) &= \frac{\partial f}{\partial q} \Big\vert_{q_0 u_0} \\ B_c(q_0, u_0) &= \frac{\partial f}{\partial u}\Big\vert_{q_0 u_0} \\ g(q_0, u_0) &= f(q_0, u_0) \end{array} \]

運動方程式

一般化座標 $ q $を使って、カートの駆動力 $ u $, カート質量 $ M $, 振り子長 $ l $, 振り子先端質量 $ m $ について運動方程式を立てます。

\[ q = [x, \theta, \dot{x}, \dot{\theta}] \]

あとの導出方法等は他所様の記事を引用して終わります

(7節まで、8節はスキップ. あと回転座標のとり方が若干違う)

Sympyモデルの準備

さて、Qiitaの記事ではラグランジアンを使っていましたが、別ルートで導出したので、こちらになります。

import sympy as sy
class cart_pole():
    def gen_rhe_sympy(self):
        g = sy.symbols('g')
        l = sy.symbols('l')
        M = sy.symbols('M')
        m = sy.symbols('m')

        q  = sy.symbols('q:{0}'.format(4))
        qd = q[2:4]
        u  = sy.symbols('u')
        
        I = sy.Matrix([[1, 0, 0, 0], 
                      [0, 1, 0, 0], 
                      [0, 0, M + m, l*m*sy.cos(q[1])], 
                      [0, 0, l*m*sy.cos(q[1]), l**2*m]])
        f = sy.Matrix([
                      qd[0], 
                      qd[1],
                      l*m*sy.sin(q[1])*qd[1]**2 + u,
                      -g*l*m*sy.sin(q[1])])
        return sy.simplify(I.inv()*f)
    
    def gen_lmodel(self):
        mat = self.gen_rhe_sympy()
        q = sy.symbols('q:{0}'.format(4))
        u = sy.symbols('u')
        
        A = mat.jacobian(q)
        #B = mat.jacobian(u)
        B = mat.diff(u)
        
        return A,B

Python

Pythonの場合はlamdifyを使って

c = CartPole()
q = sy.symbols('q:{0}'.format(4))
u = sy.symbols('u')
# calc_rhe(q, u) -> np.array
calc_rhe = sy.lambdify([q,u], c.gen_rhe_sympy())

などとしたほうが自然ですが、numbaに放り込みたいなどの事情があることもあるでしょう。そこでリファレンス通りにやってみます

from sympy.printing.pycode import pycode
c = CartPole()
pycode(c.gen_rhe_sympy())

すると、あんまりありがたくない形の出力が出ます。

'ImmutableDenseMatrix([[q2], [q3], [((1/2)*g*m*math.sin(2*q1) + l*m*q3**2*math.sin(q1) + u)/(M + m*math.sin(q1)**2)], [-(g*(M + m)*math.sin(q1) + (l*m*q3**2*math.sin(q1) + u)*math.cos(q1))/(l*(M + m*math.sin(q1)**2))]])'
  • ImmutableDenseMatrixというSympyの型そのままの出力
  • sin/cosがmath? numpyがいいな!
  • q2, q3じゃなくてq[2], q[3]と出して欲しい


さらに、リファレンスには使い方が詳しくない裏メニューのNumPyPrinterを使うとこうなります

from sympy.printing.pycode import NumPyPrinter
c = cart_pole()
NumPyPrinter().doprint(c.gen_rhe_sympy())

出力

'numpy.array([[q2], [q3], [((1/2)*g*m*numpy.sin(2*q1) + l*m*q3**2*numpy.sin(q1) + u)/(M + m*numpy.sin(q1)**2)], [-(g*(M + m)*numpy.sin(q1) + (l*m*q3**2*numpy.sin(q1) + u)*numpy.cos(q1))/(l*(M + m*numpy.sin(q1)**2))]])'

かなり良くなりましたが、

  • q2, q3じゃなくてq[2], q[3]と出して欲しい

が残ります。掘ってみるとSympyの立式で

q = sy.symbols('q:{0}'.format(4))

と定義すると、q[0] が q0 というsympyシンボルとして定義されているところから来ている様子です。

Matrixやindexingを駆使するとできるようなことも書いてある気がしないでもないですが、面倒なのでこうしました。

from sympy.printing.pycode import NumPyPrinter
c = cart_pole()
class NumPyPrinterR(NumPyPrinter):
  def _print_Symbol(self, expr):
      name = super(NumPyPrinter, self)._print_Symbol(expr)

      # 変数名が他とかぶらないのでマッチングせずに頭文字だけで判断
      if name[0] == 'q':
          name = name[0] + '[' + name[1] + ']'
      elif name[0] == 'u':
          name = 'u[0]'
      return name
NumPyPrinterR().doprint(np.squeeze(c.gen_rhe_sympy()))

これでめでたし。

C++11

Pythonで検証したあとは同じようにするだけです。C++ではEigenを使いたいのとImmutableDenseMatrixの処理が定義されていないのでもう一捻りします。

別途

Eigen::Matrix<double, 4, 1> g;
Eigen::Matrix<double, 4, 4> A;
Eigen::Matrix<double, 4, 1> B;

と定義されていることとしましょう

c = cart_pole()
class CXX11CodePrinterR(CXX11CodePrinter):
    def _print_Symbol(self, expr):
        name = super(CXX11CodePrinter, self)._print_Symbol(expr)
        
        if name[0] == 'q':
            name = name[0] + '(' + name[1] + ', 0)'
        elif name[0] == 'u':
            name = 'u(0, 0)'
        return name

print('g << ')
for expr in np.squeeze(c.gen_rhe_sympy()).tolist():
    print(CXX11CodePrinterR().doprint(expr), end=',\n')

A, B = np.squeeze(c.gen_lmodel())

print('\nA << ', end='')
for expr in A:
    print(CXX11CodePrinterR().doprint(expr), end=',\n')

print('\nB << ', end='')
for expr in B:
    print(CXX11CodePrinterR().doprint(expr), end=',\n')

あとはセミコロンとか処理してお終い。

その他

Sympyには他にも

  • julia
  • Rust

などのプリンタも用意されているそうです!

SOCPのソルバーでQPを解く話

PythonRoboticsの一例でQPのソルバーとしてECOSが使われていました。ECOSは二次錐計画問題 (SOCP)のソルバーです。 QP ⊂ QCQP(Quadratically constrained quadratic program)⊂ SOCP であって、
制約なしのQP形式の問題

\[ \begin{array}{ll} \min_{x} & \frac{1}{2} x^T Q x + q^T x + r \\ \end{array} \]

は、変数tを追加して

\[ \begin{array}{ll} \min_{x,t} & t \\ {\rm s.t.} & \frac{1}{2} x^T Q x + q^T x + r \leq t \\ \end{array} \]

と変換すればSOCPとして解ける(制約ありも同様)とのことでしたので、ちょっと実装してみます。

まず改めて、$ x \in \mathbb{R}^n $ のとき、

\[ \begin{aligned} y &\equiv \left[ x, t\right] \\ c &\equiv [0_n,1] \\ Q &= U^T U \\ Q_s &\equiv \begin{bmatrix} U & 0 \end{bmatrix} \\ q_s &\equiv [q,0] \\ \end{aligned} \]

とすると表現がこう変わります

\[ \begin{array}{ll} \min_{y} & c^T \cdot y \\ {\rm s.t.} & \frac{1}{2} y^T Q_s^T Q_s y + q_s^T y + r \leq c^T \cdot y \\ \end{array} \]

制約式をノルムの形式に変換すると

\[ \left\| \begin{matrix} \left( 1 + (q_s+c)^{T}y+r \right) / \sqrt{2} \\ Q_s y \end{matrix} \right\|_{2} \leq \frac{1-(q_s+c)^{T}y-r}{\sqrt{2}} \]

最後にCVXOPTなどでも使われるスラック変数を用いた表現に変形します。

\[ \begin{bmatrix} - (q_s + c)/\sqrt{2} \\ (q_s + c)/\sqrt{2} \\ Q_s \end{bmatrix} \cdot y + \begin{bmatrix} (1-r)/\sqrt{2} \\ (1+r)/\sqrt{2} \\ 0 \end{bmatrix} = \left[ \begin{array}{c} s_0 \\ s_1 \\s_2 \end{array} \right], s_0 \geq \left\Vert \begin{array}{c} s_1 \\ s_2 \end{array} \right\Vert_2, s_0, s_1 \in \mathbb{R}, s_2 \in \mathbb{R}^{n} \]

続けて $ G \cdot x \preceq_K h $ という形式にあわせて、

\[ \begin{bmatrix} (q_s + c)/\sqrt{2} \\ - (q_s + c)/\sqrt{2} \\ - Q_s \end{bmatrix} \cdot y + \begin{bmatrix} s_0 \\ s_1 \\ s_2 \end{bmatrix} = \begin{bmatrix} (1-r)/\sqrt{2} \\ (1+r)/\sqrt{2} \\ 0 \end{bmatrix} \]

QPをSOCPとして解く

くどいですが、以上をまとめると

\[ \begin{array}{ll} \min_{x} & \frac{1}{2} x^T Q x + q^T x + r \\ where & x \in \mathbb{R}^n \end{array} \]
というQPは、
\[ \begin{array}{ll} \min_{y} & c^T \cdot y \\ \text{s.t.} & G \cdot y \preceq_{\mathcal{K}} h \quad \left( \Leftrightarrow G \cdot y + s = h, s \in \mathcal{K} \right) \\ \text{where} & y \equiv [x\ t], t \in \mathbb{R} \\ & c = [0_{1 \times n}\ 1] \\ & Q = U^T U \\ & \mathcal{K} \equiv \{ (u_0, u_1) \in \mathbb{R} \times \mathbb{R}^{n+1} \ |\ u_0 \geq \Vert u_1 \Vert_2 \} \\ & G = \left[ \begin{array}{cc} q/\sqrt{2} & 1/\sqrt{2} \\ -q/\sqrt{2} & -1/\sqrt{2} \\ U & 0 \end{array} \right], h = \left[ \begin{array}{cc} (1-r)/\sqrt{2} \\ (1+r)/\sqrt{2} \\ 0_{n \times 1} \end{array} \right] \\ \end{array} \]
というSOCPとして解くことができます。

余談

ECOSに投げるときは、Usage from C · embotech/ecos Wikiを参照しつつ

// Nは元となるQPのxの次元
idxint l = 0; // もし制約付QPを解く場合は l>0, 配列Gに対してはconeの定義にあたる係数よりも前方に配置する
idxint ncones = 1;
idxint q[] = {N+2};

とします。

補足

$ Q = U^T U $ は、コレスキー分解です。

コレスキー分解でベンチマークを取った先日のポストも広告:Pythonと関数単位での高速化(cython, numba, swig) - 夜間飛行

参考

脱PDFからのシングルHTMLのすすめ

レポートやメモの作成に定番のMicrosoft Wordや同時編集が魅力的なgoogle documentsを使ってきたが、gitで管理できるフォーマットを試しながら、gulp+pugで作るシングルページHTMLに行き着いたのでメモする。

これまで試したもの

やりたかったこと

脱ペーパー
これまでPDF化するなど印刷前提のフォーマットを貫いてきたけどペーパーレスの時代に「紙ページに収める」という縛られ方をする必要はないのではないだろうか。 挿絵の位置で改ページがあっちいったり、こっち行ったりをどうこうしたり、 横に平べったく長い表の貼り付けたり、 章タイトルの直後に改ページが来てpandocフィルタ追加したり… そういう調節とお別れしたい。
メール添付やクラウドストレージを通した配布
クラウドストレージのプレビュー機能の充実とブラウザ間の互換性向上によって、シングルHTMLファイルがPDFの代替になりそうな気がした。
インタラクティブなコンテンツ
いったん脱PDFし、HTMLをフォーマットとして選ぶとJSとかSVGとかの採用が見えてくる。

環境の準備

Webpackも考えたけどentry pointがjsファイルになる点が今ひとつコンセプトに合わなくてgulpにした。

npm install --global gulp-cli
npm init -y
npm i -D gulp pug gulp-pug
npm i -D gulp-inline-source gulp-inline-images

構成

root
├── _template
├── dest <-- ここに出来上がったレポートが出力される
├── gulpfile.js
├── package-lock.json
├── package.json
├── img <-- レポートに使う画像を放り込んでおく
└── src <-- ここにレポートの素を置く

gulpfile

const gulp = require('gulp')
const pug = require('gulp-pug')
const inlineSource = require('gulp-inline-source')
const inlineImages = require('gulp-inline-images')

gulp.task('pug', () => {
  return gulp.src(['./src/*.pug', '!./_**/*', '!./node_modules/**/*'])
    .pipe(pug({ pretty: true }))
    .pipe(inlineSource())
    .pipe(inlineImages({ basedir: './img/' }))
    .pipe(gulp.dest('./dest/'))
})
gulp.task('watch', () => {
  return gulp.watch(['./src/*.pug', './_template/*.pug'], gulp.series(['pug']))
})
gulp.task('default', () =>  {
  return gulp.parallel(['pug'])
})

あとはgulpでもgulp watchでも…. 画像, CSSそれにjavascriptはinline化してhtmlに埋め込んでしまいます。

レポートのテンプレート

コンセプトにも関わらず、紙に印刷したいとなったらpaper.css、プレゼンしたいならreveal.jsと思いつつ、Bootstrapをcssのみで使う選択肢で作ってみた。

ライブラリを追加する場合inlineのアノテーションを忘れないように!!

doctype html
html(lang="ja")
  block config
  head
    meta(charset="UTF-8")
    meta(name="viewport", content="width=device-width, initial-scale=1")
    link(rel="stylesheet", href="https://stackpath.bootstrapcdn.com/bootstrap/4.5.0/css/bootstrap.min.css", inline)
    style.
      body {
        font-family: "Helvetica Neue", 
              "Hiragino Kaku Gothic ProN", Meiryo, sans-serif;
      }
    block css
    block scripts

    title #{title} 
body
  .container
    block content

本文の書き方

カラム数などによっては40emで折り返すようにしておくと読みやすいかも。

extends ../_template/template.pug

block config
  - var title='私のレポート'

block css
  style.
    .alert-trim {
      display: inline-block;
    }
    .col p {
      width: 40em;
    }

block content
  .row
    .col
      h1 報告
      h2 はじめに
      p 今月は…

Boilerplate

ボイラープレート用意した

Chachay/pug_boilerplate: My Boilerplate of Single HTML reports