coord_desc_NNLS.jl

function NNLS_CD(X, Y; ϵ = 1e-9, max_iters = 300)
    rows,vars = size(X)
    XTX,XTY = transpose(X) * X, transpose(X) * Y
    x  = zeros(vars)
    μ = -XTY
    Hxf = similar(XTY)
    @inbounds for iter in 1:max_iters
        Hxf = XTX * x - XTY
        all(>=(-ϵ), Hxf) && break
        for v in 1:vars
            initial = x[v]
            x[v,1]  = max(x[v] - μ[v] / XTX[v,v], 0.0)
            ∇       = x[v] - initial
            μ      .+= ∇ .* @views XTX[:,v]
        end
    end
    
    return x
end

variables = 100
X    = randn(200, variables)
beta = rand(variables, 1)
Y    = X * beta
r = NNLS_CD(X, Y; ϵ=1e-9, max_iters=300)
sum(abs,beta .- r)

using BenchmarkTools
#@simd - 6.459ms
#no simd - 6.420ms
#faster check: 6.339
#@views: 4.530 ms
#better?: 3.576
#wow?  2.031 ms 
@btime NNLS_CD(X, Y; ϵ=1e-9, max_iters=300)