2015年5月22日金曜日

Julia の発展途上なところ

この前の簡単な内点法などをしてみて、以下のところが数理最適化に使うには Julia はまだまだ発展途上な印象を受ける。

(1) 遅い
C, C++ などと比較しても遜色ない速さ、とうたっているが、これはある意味で当然で、そもそも実行前にコンパイルしている(インタープリターではない)。

また、実際にプログラミングをしてみると
using PyPlot
など頻出する using の遅さが地味に結構痛い。

ちなみに、これについては、PyPlot などを ~/.juliarc でロードしておけば、それなりに快適にできる。Matlab が起動に時間がかかっている作業を Julia でも行っているようなもので、起動に時間をかけることで、そのあとの作業を効率化できる。


(2) Profiling がカウンターであって秒で計測できない。
たとえば 1/1000 秒ごと、と指定してもサンプルカウントをしているので、かなり荒い情報しか得られない。だいたいのボトルネックを把握するには有効であるが、大まかなチューニングしか行えない。
しかも @time などの時間測定のコマンドを組み入れた関数を @profile で測定すると、サンプルカウントがうまくいかないことがある。

(3) IDE の機能がいまひとつ。
Juno を使ってみたが、数理最適化に多い反復計算には向いていなかった。現状では、emacs とコマンドラインが、一番便利。
ただし、Julia の REPL は結構強力であるので、コマンドラインからの利用でもそれなりに行けたりする。



ここに挙げたうちで一番重要そうなのは (3) で、たぶん IDE に Python の spyder なみの環境ができれば Julia の利用可能性はぐっと高まるかと思う。



2015年5月21日木曜日

TeXインストーラ3での「mktexlsr の実行に失敗しました」を修正する

TeX インストーラ3で TeX を一括で更新した時に、なぜか
\usepackage{color}
などで
! LaTeX Error: File `color.sty' not found.
と怒られるときがある。
これは、TeX インストーラ3 の更新の時に mktexlsr の実行に失敗しているためのようで、abtexinst_log.txt にもエラーメッセージとして記録されている。

これを修正するには、コマンドプロンプトを「管理者権限」で開いて
$ mktexlsr
と実行する。
(「管理者権限」でなくても mktexlsr は実行できて成功するが、ファイルには保存されないので実質的に「管理者権限」が必要である。)


2015年5月14日木曜日

実用的な最適化とは?

シミュレーテッドアニーリングで max cut を解く、という論文が arXiv の

Solving maximum cut problems by simulated annealing

で読むことができる。
シンプルな実装(ソースコードも途中に書いてあるぐらい)になっているが、数値結果を見てみると、それなりに良い結果になっているようにも見える。


最適化の論文には難しい理論が書かれているマニアックなものが多くて、そういうのが多いのも面白いが、こういう simple is best みたいなことも研究の方向としてはあるべき一つなのかと思ったりもする。
この論文がどこの雑誌に採択されるか、気になるところだ。

2015年5月13日水曜日

Julia で戸惑うところ (Julia と Matlab の違い)

Julia は Matlab に似たように書けるので、ついつい同じように書いてハマることがあるので、そういったことについてのメモをしておく。

(1) 変数の deepcopy に気をつけろ!
たとえば、Matlab で
x = [3 4 5];
y = x;
y(2) = 7;
としても、x は変わらず x = [3 4 5];である。
しかし、Julia で
x = [3 4 5];
y = x;
y[2] = 7;
とすると、y につられて x も x = [3 7 5] に変化する。
要するに、Matlab では y = x で別の配列が準備されて配列の中身もコピーされるが、Julia の y = x は単純に「 x という配列に y という名前でもアクセスできる」ということで中身も同じである。
Matlab のようにするには、
y = deepcopy(x)
とする。


(2) integer と float に気をつけろ!
Matlab の場合は基本的に float (double) で計算されているが、Julia の場合は厳密に違う型として認識されている。
たとえば、
x = [3 4 5]
とすると出来上がる型は 1x3 Array{Int64,2}: である。
したがって、
x[1] += 2.3
とすると、float の値となる 3 + 2.3 = 5.3 を x[1] に代入しようとして
ERROR: InexactError()
で失敗する。

このときには、
x = float([3 4 5])
と float をかませておけば回避できる。


(3) max と maximum の違いに気をつけろ!
Matlab の max 関数は、Julia では max と maximum に機能の違いで分割されている。
max は要素を引数として並べることできる状態、maximum なら配列が引数となる。
たとえば、

max(2,3) と直接 2と3が引数となる場合は max, maximum([2,3]) と配列[2,3]が引数ならmaximum である。
逆に max([2,3]), maximum(2,3) とするとエラーで失敗する。




2015年5月11日月曜日

Julia の面白いところ(数理最適化のプログラミングから見て)

この前、内点法を実装してみて分かった、面白いところを列挙してみる。

1. UTF-8 の文字列を変数に使える。例えば
  α = β + γ
なんて式が可能。

2. A ? B : C が使える。(C っぽい)

3. リスト内包が使える。(python っぽい)
これは、SDP のときに便利で、例えば内積の A_k \bullet X というのを A_1,...,A_m で行いたければ
float([ sum(A[:,:,k].*X) for k=1:m ])
とすれば一行で作れる。
(float をかけているのは、リスト内包で上のリストを作るとなぜか Any のリストになってしまって、その後に足し算引き算などができなくなってしまうので、float に型を強制変更している。)

4. Cell にできる。(Matlab っぽい)
たとえば、
Z = cell(3,2)
Z[1,2] = [3 4; 5 6]
Z[3,1] = spzeros(7,8)
のように違うデータを入れ込むことができる。


3 と 4 は SDP のときには相当便利で、これがあれば SDPA-M 入力と SeDuMi 入力の両方のアドバンテージを一つのデータでこなすことができる。

2015年5月5日火曜日

Julia で SDP の内点法を実装してみた

次のようなソースファイルを sdpa01.jl として保存しておいて
julia > include("sdpa01.jl")
とすると内点法で SDP を解ける。例題は、SDPA にある example1.dat-s を解いている。


# input data
m = 3;
n = 2;
A = zeros(n,n,m)
A[:,:,1] = [10 4; 4 0]
A[:,:,2] = [0 0; 0 -8]
A[:,:,3] = [0 -8; -8 -2]
b = [48, -8, 20]
C = [11 0; 0 -23]

# parameters
ε = 1.0e-6
γ = 0.9
β = 0.1
λ = 1.0e+4
maxIter = 30

# initialization
iter = 0
X = λ*eye(n);
Y = λ*eye(n);
z = zeros(m,1);

M = zeros(m,m); # Schur complement
# main loop
while iter < maxIter
iter += 1
# compute residual
p = b - float([sum(A[:,:,k].*X) for k=1:m])
   # float is necessary here
D = C - Y - sum([ A[:,:,k]*z[k] for k=1:m])
max_p = maximum(abs(p));
max_D = maximum(abs(D))
μ = sum(X.*Y)/n
residuals = max(max_p, max_D, μ);
if  residuals <= ε
@printf "iter = %3d, r_p = %.2e, r_d = %.2e, μ = %.2e, α = %.2e\n" iter max_p max_D μ α
println("Converge !!")
break
end

Yinv = inv(Y)
# build Schur complement
for j=1:m
XAY = X*A[:,:,j]*Yinv
M[:,j] = float([sum(A[:,:,i].*XAY) for i=1:m])
end

# search direction
dz = M \ (p - float([ sum(A[:,:,k].* (β*μ*Yinv - X - X*D*Yinv)) for k=1:m]) )
dY = D - sum([ A[:,:,k]*dz[k] for k=1:m ]);
dX = β*μ*Yinv - X - X*dY*Yinv
dX = (dX + dX')/2;

# step length
Lx = chol(X,:L); Ly = chol(Y,:L)

(Dx, Vx) = eig(inv(Lx)*dX*inv(Lx)')
eigx = minimum(Dx)
αp = eigx < 0 ? -1.0/eigx : 100.0
(Dy, Vy) = eig(inv(Ly)*dY*inv(Ly)')
eigy = minimum(Dy)
αd = eigy < 0 ? -1.0/eigy : 100.0

α = min(γ*αp, γ*αd, 1.0)
# update
X += α*dX
Y += α*dY
z += α*dz
@printf "iter = %3d, r_p = %.2e, r_d = %.2e, μ = %.2e, α = %.2e\n" iter max_p max_D μ α
end

if iter <= maxIter
println("--Optimal Solution--")
println("X = ", X)
println("Y = ", Y)
println("z = ", z)
else
println("Failed to converge during ", iter, " iterations.")
end