Hi,
I’m writting a simple Dykstra sheme for projecting a vector start
onto equality constraints Ax = b
and inequelity constraints x >= 0
. It works great, but I wander if it could be faster ?
Here’s the current version :
using Random
Random.seed!(11)
N,M = 30,50
A = big.(rand(N,M));
x0 = big.(rand(M));
b = A*x0;
# check :
A*x0 == b;
start = big.(randn(M));
get_AiAA(A) = A'inv(A*A')
proj_eq(x,A,b,AiAA) = x .- AiAA*(A*x-b)
proj_ineq(x::T) where T = x > T(0) ? x : T(0)
function naive_version(start,A,b)
AiAA = get_AiAA(A)
x = deepcopy(start)
old_x = zero(x)
y = zero(x)
p = zero(x)
q = zero(x)
# iter = 0
while !(old_x ≈ x)
old_x .= x
y .= proj_eq(x.+p,A,b,AiAA)
p .= x .+ p .- y
x .= proj_ineq.(y .+ q)
q .= y + q - x
# iter += 1
end
# println("Converged in $iter iterations")
return x
end
function update!(old_x,x,p,q,ProjP,Projq)
old_x .= x
x .= x .+ p
p .= ProjP*x .- Projq
x .= x .- p .+ q
q .= x
x[x.<0] .= 0
q[q.>0] .= 0
end
function Dykstra(start,A,b)
# Deal with the linear projection
AiAA = A'inv(A*A')
ProjP = AiAA*A
Projq = AiAA*b
# Dykstra :
x = deepcopy(start)
old_x, p, q = zero(x), zero(x), zero(x)
while !(old_x ≈ x)
update!(old_x, x, p, q, ProjP, Projq)
end
return x
end
naive_solution = naive_version(start,A,b);
solution = Dykstra(start,A,b);
using Test
using BenchmarkTools
@test solution ≈ naive_solution
@test all(solution .>=0)
@test A*solution ≈ b
@btime Dykstra($start,$A,$b);
I included a naive version that is known to be the right algorithm to compare to. On my machine it takes 2.1s
for these array sizes, but it gets slow for bigger arrays.
It seems to be type-stable. I have concerns regarding broadcasting : could’nt it be faster with loops in the update!
function ?