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 ?