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