数値実験の続き

CPLEX で混合整数計画問題を解き始めて、はや１週間。
結局、CPUを 1500% 使用し続けたが上界と下界のギャップが8%程度までしか接近せず、n = 1000 の問題は断念。

もっと小さい問題、ということで n = 200 にして再度挑戦中。10% 程度までは比較的簡単だが、やはり本当の勝負はそこからの様子。

cplex が思っていた以上に時間がかかる

SOCP 付きの整数計画問題で厳密解がどうなっているかを知るために、あまり深く考えずに cplex にお任せしたら、想定以上に時間がかかっていることに気がついた。

変数の数は、0-1 変数が 1000, 連続変数も 1000 で、0-1 変数を 0<=x<=1 に緩和した SOCP は 0.1 秒もかからずに解けるので、それほど大きい問題ではないかと思い、そのまま cplex に丸投げしてある。

今のところ CPU 使用率が 1500% 程度を 24 時間ぐらいやってみて、上界と下界の差が 0.2% から縮まらない状態である。

分枝限定法での限定が上手くできないタイプの問題なのかもしれないが、もう少し定式化もいじってみようと思う。

Ipopt の Matlab インターフェースのコンパイルで No such file or directly と言われる

Ipopt を試すためにコンパイルしてみたが Matlab インターフェイスでコンパイルがうまく行かなかった。

より具体的には、
CoinIpopt/Ipopt/contrib/MatlabInterface/src
のディレクトリで
$ make
として ipopt.mexa64 を生成しようとすると、以下のエラーで停止する。

g++: error: matlabexception.o: No such file or directory
g++: error: matlabfunctionhandle.o: No such file or directory
g++: error: matlabjournal.o: No such file or directory
g++: error: iterate.o: No such file or directory
g++: error: ipoptoptions.o: No such file or directory
g++: error: options.o: No such file or directory
g++: error: sparsematrix.o: No such file or directory
g++: error: callbackfunctions.o: No such file or directory
g++: error: matlabinfo.o: No such file or directory
g++: error: matlabprogram.o: No such file or directory
g++: error: ipopt.o: No such file or directory

make のときのエラーを確認していったところ、.libs サブディレクトリにこれらのファイルがあることがわかったので、以下のように対応してみた。

$ cd CoinIpopt/Ipopt/contrib/MatlabInterface/src

$ cp .libs/*.o .

$ make

これで無事に ipopt.mexa64 が生成できる。

IguanaTeX で日本語が通るようになった！

http://karat5i.blogspot.jp/2015/04/iguanatex-powerpoint-tex.html
に書いてある通り、 IguanaTeX で日本語が通るようになった。
自分のところでも試してみたが、正常に動作している。

今までは日本語をうまく使えていなかったが、これで platex でコンパイルできるものをPowerPoint にも使えるようになり、とても便利になった。

素晴らしい。

学術雑誌で投稿ができない状態になっているものがある

数理最適化関係の学術雑誌で、投稿数が増えすぎて新規投稿を reject している雑誌がある。
別の雑誌に投稿するように、という意味での reject になっている。
ただ、ホームページにそういったことが書かれてはいないので、投稿するまではその雑誌がそういう状態になっていることを把握するのは難しそうだと思う。

それにしても、査読などをたまに引き受けたりもするが、やっぱり査読者を見つけてくるという Editor や Associate Editor のお仕事は相当大変なものだろうから、そのあたりはほんと感謝だな、と思ったりする。

Julia でループ展開しても、あまり速くない？

疎行列と密ベクトルの積を計算するような以下のコードを書いてみた。

==== Ax.jl =====

function Ax(n, density, seed)
    srand(seed)
    println("Generating A")
    time_A = @elapsed @time A = sprand(n,n,density)
    x = rand(n,1)
    println("Computing b = A*x")
    time_b = @elapsed @time b = A*x

    println("Computing b2 = A*x")
    time_b2 = @elapsed @time begin
        b2 = zeros(n,1)

        for j=1:n
            const xj = x[j]
            for index_i = A.colptr[j]:(A.colptr[j+1]-1)
                b2[A.rowval[index_i]] += A.nzval[index_i]*xj
            end
        end
    end

    println("norm(b-b2) = ", norm(b-b2))

    println("time (A) = ", time_A, " time (b) = ", time_b, " time(b2) = ", time_b2)

    return A, x, b, b2

end

==============

実行結果は以下の通り。

julia> include("Ax.jl"); (A, x, b,b2) = Ax(2000,0.8,1024); (A,x,b,b2) = Ax(2000,0.8,1024);
Generating A
elapsed time: 0.173287485 seconds (76853144 bytes allocated, 24.33% gc time)
Computing b = A*x
elapsed time: 0.011240257 seconds (16064 bytes allocated)
Computing b2 = A*x
elapsed time: 2.677362552 seconds (498951072 bytes allocated, 9.45% gc time)
norm(b-b2) = 0.0
time (A) = 0.173453157 time (b) = 0.011356873 time(b2) = 2.677506154
Generating A
elapsed time: 0.141352727 seconds (76853144 bytes allocated, 12.73% gc time)
Computing b = A*x
elapsed time: 0.010914271 seconds (16064 bytes allocated)
Computing b2 = A*x
elapsed time: 2.708701738 seconds (498951072 bytes allocated, 9.03% gc time)
norm(b-b2) = 0.0
time (A) = 0.141499582 time (b) = 0.011023018 time(b2) = 2.708858507

b の計算はもともとあるもので、b2 のほうがループを展開して計算している。

ちなみに、２回実行しているのは、Julia のコンパイルで時間がかかると遅くなるかもしれない、ということで２回実行してみた。

うーむ、やはり numpy のように、できるだけ行列やベクトルの演算に持ち込んだ方が

速いのだろうか？

そうそう、@elapsed, @time マクロは begin - end で挟むと、その間の時間を測定してくれる。

@elapsed マクロのほうは返り値で経過時間を返してくれるのでそれを代入しているが、@time マクロのほうが細かい情報が表示される。

Julia で簡単な Newton 法の実装

関数 f_1(x) = x[1]^4 - 2x[1]x[2]+ 2x[2]^2 - 4x[1] + 2x[3] についての Newton 法を実装してみると以下のようになり短いコードで書ける。

面白い点としては、

1. C や Matlab なら 2*x[1]*x[2] と書かなければいけないコードが 2x[1]x[2] と * を省略して書ける
2. UTF-8 でファイルを作っていれば、∇も普通の文字として扱える (全角空白が紛れ込むと、かなり厄介だけど)

実行するには
julia> include("newton.jl")
とする。
ちなみに、このコードなら短いので 15 分程度で書ける。

==== newton.jl ====

# -*- coding:utf-8 mode:Julia -*-

function f_1(x)
    return x[1]^4 - 2x[1]x[2]+ 2x[2]^2 - 4x[1] + 2x[3]
end
function ∇f_1(x)
  return [4x[1]^3-2x[2]-4; -2x[1]+4x[2]+2]
#    return [4*x[1]^3-2*x[2]-4; -2*x[1]+4*x[2]+2]
end
function ∇2f_1(x)
    return [12x[1]^2 -2; -2 4]
end


ε = 1.0e-8
k = 0
x = [1; 1]
norm_∇f = 0.0
while k < 100
    ∇f = ∇f_1(x)
    norm_∇f = norm(∇f)
    if norm_∇f < ε
        println("Convergence!!")
        break
    end
    ∇2f = ∇2f_1(x)
    @printf "k = %d : x = [%.2e, %.2e], norm(∇f) = %.2e\n" k x[1] x[2] norm_∇f
    x = x - ∇2f \ ∇f
    k = k+1
end
@printf "k = %d : x = [%.2e, %.2e], norm(∇f) = %.2e\n" k x[1] x[2] norm_∇f

==========