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

==========

内点法で新しい結果

arXiv をチェックしていたら、内点法の新しい結果が出ていて、なんと infeasible interior-point method のオーダーが低くなっている。

内点法や線形計画問題の理論的な性質は 2000 年頃までの研究でほとんどなされており、それ以降はソフトウェア開発などが中心で理論的な前進はあまりない状況だった。
理論的な論文に関しても、「今までと同じオーダーの解析ができました」が主流であったが、今回の論文は違うようだ。

非常に面白そうなので、あとで読んでみることにしよう。

Julia で最小自乗法を組んでみる

簡単な最小自乗法なら、以下のような感じでプログラムを書ける。
Matlab にそっくりな感じ。

==== least1.jl =====
function leastsquare(m, n, seed)

    if m < n
        error("Too small m")
    end
    srand(seed)
    A = rand(m,n)
    b = rand(m)

    xhat = (A'*A) \ (A'*b)
    e = A*xhat - b
    norm_e = norm(e)
    @show A
    @show b
    @show e
    @printf "norm_e = %.3e\n" norm_e

end
===============

実行するには、
julia> include "least1.jl"
julia > leastsquare(5, 2, 1024)
といった感じ。

なんだか良く分からないバグが取れた

Ｃ＋＋で組んでいたときに、

if (data[419] ! = data2[419]) {
   printf("diff = %d\n", data[419]-data2[419]);
}

としたら、なぜか diff = 0 と表示される、というバグが発生。
なぜなのか良く分からないが、data, data2 にデータを入れる前に全部の要素を0で初期化するとバグが取れた。


おそらくメモリリークしているのかと推測されるが、他のところがあまりに単純なコード過ぎてメモリリークを発見できず。
しかも、Linux だと発生せずに Windows 上だと発生する、というたちの悪いバグになっており、valgrind などでもうまく特定できず。

なんかよく分からないけど、とりあえずバグが取れたのでいいにしておこうかと思う。

TeX でå (a の上に丸がついたもの)を使うマクロ

TeX で a の上に丸がついた å を使うとき、本来であれば、
\r{a}
とすればＯＫだが、\r を別のマクロで使っていたりすると、これが困る。

この回避方法は、automake の texinfo.tex を見ていて分かったが、
\def\ringaccent#1{{\accent23 #1}}
とマクロを定義してから
\ringaccent{a}
とすると回避できる。


ドイツ語のウムラウトなども\accent24 などで定義されていたりするので、これをうまく使うと \o などを別のマクロに使ったりもできて便利でもある。

Julia で Rosenbrock 関数の最急降下法

練習用に実装をしてみた。
以下のファイルを rosen.jl として保存して、
include("rosen.jl"); test01(3);
とすれば、Rosenbrock の n=3 のときを最小化してくれる。
（ただし、4154 反復かかる。）

やはり、Julia で便利なのは
1. リスト内包ができる
2. UTF-8 の文字が変数名に利用できる（∇も利用できる）
3. @show マクロが使える
のあたりかと思う。
Rosenbrock の関数式を調べるところからはじめて、この文章を書くまでで、だいたい1時間といったところ。

そういえば、C 言語なら x == 0 || y == 0 のようなところの || の前に改行が入ってもOKだが、Julia の場合は NG になる様子。

== ここから rosen.jl ==

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

function Rosen(x)
    sum([ (1-x[i])^2 + 100*(x[i+1] - x[i]^2)^2 for i = 1:length(x)-1])
end

function ∇Rosen(x)
    N = length(x)
    [[-2*(1-x[i])+2*100*(x[i+1]-x[i]^2)*(-2*x[i]) for i=1:length(x)-1]; 0] + [0; 2*100*[x[i+1] - x[i]^2 for i=1:length(x)-1]]
end

function minimizer(n, f, ∇f)
    count_f   = 0
    count_∇f = 0
    η = 0.05

    x = zeros(n,1)
    current_f = f(x)
    count_f += 1
    max_iter = 10000
    iter = 1
    while iter <= max_iter
        gf = ∇f(x)
        norm_gf = norm(gf)
        count_∇f += 1
        α = 1.0
        next_f = 0.0
        next_x = 0
        while true
            next_x = x + α*(-gf/norm_gf)
            next_f = f(next_x)
            count_f += 1
            if next_f <= current_f - η*α*norm_gf || α < 1.0e-10
                break
            end
            α *= 0.5
        end
        @printf("iter: %d, f: %.4e, gf = %.2e, α = %.2e\n", iter, current_f, n\
orm_gf, α)
        if norm_gf < 1.0e-4
            println("Norm of gradient is small")

            break

        end

        if α < 1.0e-10

            println("Step length is too short")

            break

        end

        if abs(current_f - next_f)/max(abs(current_f),1.0) < 1.0e-10

            println("No improvement of f")

            break

        end

        current_f = next_f

        x = next_x

        iter += 1

    end

    x

end

function test01(n)

    minimizer(n, Rosen, ∇Rosen)

end