# Speedup of \`\\\` solver bottleneck

**URL:** <https://discourse.julialang.org/t/speedup-of-solver-bottleneck/11279>\
**Category:** Performance\
**Tags:** question, indexing, linearalgebra\
**Created:** [May 31, 2018, 5:47am UTC](https://discourse.julialang.org/t/speedup-of-solver-bottleneck/11279 "2018-05-31T05:47:19Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![stst](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stst/32/17252_2.png) [@stst](https://discourse.julialang.org/u/stst)\
**Post date:** [May 31, 2018, 5:47am UTC](https://discourse.julialang.org/t/speedup-of-solver-bottleneck/11279/1 "2018-05-31T05:47:19Z")

</div>

Below is my implementation of the block principal pivoting algorithm for solving a non-negative linear least squares problem. This function is the bottleneck in my code (running on v0.6.2), which calls this in an inner loop with square matrices `AtA` of size 200-600 and appropriate vectors `Atb` and `x0` as inputs.

Is there any way to speed this up?

Based on profiling a Matlab version of this, it appears that essentially all the time is spent solving the linear system: `x[F]=AtA[F,F]\Atb[F]`, where `F` is a mask that changes significantly every iteration. I have tried prepending `@views` and/or `@inbounds`, without any effect on the timings (measured using `@btime`). Also, `@code_warntype` does not indicate any type instability.

Any advice would be much appreciated!

Here is the full function:

```julia
function nnls_bpp(AtA,Atb,x0)

    n = length(Atb)
    F = x0 .> 0

    x = zeros(n)
    x[F] = AtA[F,F]\Atb[F]
    y = AtA*x - Atb
    xFneg = (x .< 0) .& F
    yGneg = (y .< 0) .& .!F
    nInfeasible = sum(xFneg .| yGneg)

    p = 3
    t = n+1
    while nInfeasible > 0

        if nInfeasible < t
            t = nInfeasible
            p = 3
            F[xFneg] = false
            F[yGneg] = true
        elseif p >= 1
            p = p - 1
            F[xFneg] = false
            F[yGneg] = true
        else
            idx = findlast(xFneg .| yGneg)
            F[idx] = .!F[idx]
        end

        x = zeros(n)
        x[F] = AtA[F,F]\Atb[F]
        y = AtA*x - Atb
        xFneg = (x .< 0) .& F
        yGneg = (y .< 0) .& .!F
        nInfeasible = sum(xFneg .| yGneg)

    end

    return x

end

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 31, 2018, 5:52am UTC](https://discourse.julialang.org/t/speedup-of-solver-bottleneck/11279/2 "2018-05-31T05:52:36Z")

</div>

Have you tried using other solvers like IterativeSolvers.jl yet? Are there special matrix properties to exploit?

---

<div class="post-metadata">

**Author:** ![stst](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stst/32/17252_2.png) [@stst](https://discourse.julialang.org/u/stst)\
**Post date:** [May 31, 2018, 6:22am UTC](https://discourse.julialang.org/t/speedup-of-solver-bottleneck/11279/3 "2018-05-31T06:22:09Z")

</div>

IterativeSolvers.jl does not seem to support non-negativity constraints needed for this NNLS problem.

I have considered the old Lawson-Hanson active-set algorithm, as well as the somewhat more recent “Fast NNLS” by Bro-De Jong. The block prinicipal pivoting algorithm turns out to be an order of magnitude faster than these.

The matrix `AtA` is symmetric and appears to be positive definite.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 31, 2018, 6:24am UTC](https://discourse.julialang.org/t/speedup-of-solver-bottleneck/11279/4 "2018-05-31T06:24:06Z")

</div>

> [@stst](#):
>
> The matrix `AtA` is symmetric and appears to be positive definite.

If it’s SPD then use `cg` from IterativeSolvers.jl. Or you can try wrapping `AtA` in a `Symmetric` before `\` to speed some operations up.

> [@stst](#):
>
> IterativeSolvers.jl does not seem to support non-negativity constraints needed for this NNLS problem.

And neither does `\`?

---

<div class="post-metadata">

**Author:** ![stst](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stst/32/17252_2.png) [@stst](https://discourse.julialang.org/u/stst)\
**Post date:** [May 31, 2018, 6:40am UTC](https://discourse.julialang.org/t/speedup-of-solver-bottleneck/11279/5 "2018-05-31T06:40:54Z")

</div>

Thanks for the suggestions! (Sorry, I misunderstood you - thinking to use `cg` to replace my entire function.)

I just tried `cg` from `IterativeSolvers`. Unfortunately, it is 3-5 times slower than `\` for my problem.

Surprisingly, replacing `AtA[F,F]` by `Symmetric(AtA[F,F])` results in a 20% slowdown.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [May 31, 2018, 7:47am UTC](https://discourse.julialang.org/t/speedup-of-solver-bottleneck/11279/6 "2018-05-31T07:47:41Z")

</div>

> [@stst](#):
>
> IterativeSolvers.jl does not seem to support non-negativity constraints needed for this NNLS problem.

Try `Optim.jl` to solve the whole problem. Use [Fminbox{ConjugateGradient}](http://julianlsolvers.github.io/Optim.jl/stable/user/minimization/#box-minimization).

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [May 31, 2018, 12:29pm UTC](https://discourse.julialang.org/t/speedup-of-solver-bottleneck/11279/7 "2018-05-31T12:29:00Z")

</div>

In case you want to try Lawson-Hanson, see [https://github.com/rdeits/NNLS.jl](https://github.com/rdeits/NNLS.jl) for a good implementation (a direct port of the old Fortran code that doesn’t allocate and is fully generic).

---

<div class="post-metadata">

**Author:** ![Juser](https://avatars.discourse-cdn.com/v4/letter/j/34f0e0/32.png) [@Juser](https://discourse.julialang.org/u/Juser)\
**Post date:** [May 31, 2018, 8:51pm UTC](https://discourse.julialang.org/t/speedup-of-solver-bottleneck/11279/8 "2018-05-31T20:51:54Z")

</div>

Perhaps there is some sort of analog to FastMath for large matrix multiplication/division that solves x\y much more approximately. You probably don’t need a huge amount of precision (possibly ever, but certainly not until your final round of iterations). You might also be able to talk with BLAS directly and try to exploit information about your matrix (e.g. symmetry).

---

<div class="post-metadata">

**Author:** ![stst](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stst/32/17252_2.png) [@stst](https://discourse.julialang.org/u/stst)\
**Post date:** [June 1, 2018, 3:33am UTC](https://discourse.julialang.org/t/speedup-of-solver-bottleneck/11279/9 "2018-06-01T03:33:19Z")

</div>

Thanks for all the hints!