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

0 件のコメント:

コメントを投稿