簡単そうに見えて、実は難しい、最適化問題

最小自乗問題
$$ \min \frac{1}{2} || Ax - b ||^2 \quad \mbox{s.t.} \quad x \in R^n $$
は、A の行が一次独立なら簡単に解ける問題。

次に、
$$ \min \frac{1}{2} || x - x^0 ||^2 \quad \mbox{s.t.} \quad x \in R^n,\ x \ge 0$$

と各変数が非負に制限されるときも簡単に解けて、解としては

$$ x^* = \max\{x,0\} $$
である。

ところが、この２つが組み合わさると
$$ \min \frac{1}{2} || Ax - b ||^2 \quad \mbox{s.t.} \quad x \in R^n,\ x \ge 0$$

とたんに、簡単には解けなくなる。

http://math.uni-graz.at/mobis/publications/SFB-Report-2012-016.pdf

のように内点法ベースなら解けるのであるが、もっと簡単に解けそうにも見えて、そういった解法が見当たらない。

この「簡単そうなのに、より適した解法が案外見当たらない」というところが、面白い。

LPの今までとは違う解法

最近の論文をチェックしていたところ、LP のシンプレックス法や内点法とは異なる解法についての論文があった。
ぱらぱらっとめくったところでは、最適解が得られることなどは証明されている様子だが、どれだけの反復回数で収束するのか、のあたりはまだ検討されていない様子。ただ、「反復回数」の扱いも内点法と同じような基準でいいのかどうかもよく分からないので、単純な比較はできないかも知れない。

流れとしては、シンプレックス法よりも内点法に近い流れなので、中心パスなどとどう関係しているのか、などを考えてみると面白いかと思う。

unit simplex への射影

unit simplex、つまり
$$ \sum_{i=1}^n x_i= 1, x \ge 0 $$
への射影の方法が以下の論文に掲載されている。

A projected gradient method for optimization over density matrices
D.S. Gonçalves, M.A. Gomes-Ruggiero & C. Lavor
http://www.tandfonline.com/doi/abs/10.1080/10556788.2015.1082105?journalCode=goms20

Dykstra などの交互射影法とは異なり、有限回の計算で解が計算される点が面白い。

ところで、この方法を変形して、
$$ \sum_{i=1}^n x_i= 1 $$
にまずは射影してから計算すると効率が上がったりするのだろうか？
そのあたりを考えてみるのも面白そうではある。

著者が亡くなっても掲載される論文

Mathematical Programming Computation の Volume 7, Number 3 の Preface にいろいろと事情が掲載されている。

非常に意義深いことだな、と思う。

Matlab ソースを PDF に一括変換する

Matlab のソースを Matlab がインストールされていないパソコンで閲覧したいときには一度 PDF に変換したりすると便利だが、 Matlab の機能の publish だと HTML に変換して PDF に変換するのを一つ一つ行う必要があって、これがなかなかに面倒。

ということで、以下のようなシェルスクリプトで一括変換する。

for i in `find . | grep m$`
do
    LANG=C a2ps -1 --pro=color $i -o - | ps2pdf - $i.pdf
done

ここで、 LANG=C を入れているのは、a2ps で日本語が文字化けするので日付などが日本語だと文字化けしてしまうため。
また、--pro=color とすると、出力がカラーになるので、見やすくなる。

対角優位行列とゲルシュゴリンの定理

「対角優位行列は半正定値行列である」

という命題が成り立つ方法を考えていたが、ゲルシュゴリンの定理を使ったら、あっというまに示せることがわかった。

やっぱり、いろんな定理を知っている、ということも大事なものである。

Julia で多倍長計算

Julia でも多倍長計算できる、ということでチェックしてみた。

基本的には BigFloat を使えば簡単に行うことができて、with_bigfloat_precision を使えば特定の範囲だけ精度を変更できる。

with_bigfloat_precision(30) do
    x = BigFloat(0.1)
    i = 0
    while true
        i += 1
        x = x/2
        y = e^(-x^2)
        print(i)
        print(",")
        print(x)
        print(",")
        println(y)
        if y==1
            break
        end
    end

end

このようにすると、30ビット分の精度で計算できるし、30を1000にすれば、もちろんもっと精度が出せる。

ちなみに、

X = [BigFloat(rand()) for i=1:5, j=1:5]

とすれば、5x5 の行列なども作れる。

足し算や inv など基本的な操作はすでに行えるが、eig や chol のあたりはまだ実行できない様子である。

行列周りが完全ではないけど、いろいろとできる範囲はそれなりに広くて、なかなかに面白そうな機能だ。