# Optimzing many linear solves

**URL:** <https://discourse.julialang.org/t/optimzing-many-linear-solves/100304>\
**Category:** Performance\
**Created:** [June 13, 2023, 10:50pm UTC](https://discourse.julialang.org/t/optimzing-many-linear-solves/100304 "2023-06-13T22:50:00Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![RobertGregg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robertgregg/32/22105_2.png) [@RobertGregg](https://discourse.julialang.org/u/RobertGregg)\
**Post date:** [June 13, 2023, 10:50pm UTC](https://discourse.julialang.org/t/optimzing-many-linear-solves/100304/1 "2023-06-13T22:50:00Z")

</div>

I have a function that runs many linear solves on subsets of columns from a matrix I give it. Something like:

```julia
function score(data, parents, child)

    #Subset columns from the data matrix
    @views begin
        X = data[:, parents]
        y = data[:, child]
    end

    #Solve the system
    b = X\y

    #Get the predicted values
    ŷ = X*b

    #return mean square error
    return sum( (yᵢ - ŷᵢ)^2 for (yᵢ, ŷᵢ) in zip(y,ŷ)) / length(y)
end

```

An example of calling this function:

```julia
data = rand(15,10)
parents = [1,2,3]
child = 4
result = score(data, parents, child)

```

My question is whether there are any further optimizations I can add to make this function as performant as possible. In my case, the size of the `parents` vector can change which I think makes things tricky.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 13, 2023, 11:16pm UTC](https://discourse.julialang.org/t/optimzing-many-linear-solves/100304/2 "2023-06-13T23:16:00Z")

</div>

If you do the QR factorization of `data` (assuming that is full column rank), you can re-use it for repeated least-square solves with a subset of the columns. See e.g. [GitHub - mpf/QRupdate.jl: Column and row updates to "Q-less" QR decomposition, including stable least-squares solves](https://github.com/mpf/QRupdate.jl/tree/master#user-content-deleting-columns)

---

<div class="post-metadata">

**Author:** ![RobertGregg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robertgregg/32/22105_2.png) [@RobertGregg](https://discourse.julialang.org/u/RobertGregg)\
**Post date:** [June 14, 2023, 1:59am UTC](https://discourse.julialang.org/t/optimzing-many-linear-solves/100304/3 "2023-06-14T01:59:59Z")

</div>

Oh this looks really interesting. In my case there’s no reason why the data wouldn’t be full rank. I tried implementing this method and it gave a little bit more than a 2x speed up 👍. However, occasionally it gives different/worse solutions (relative to `X\y`). It looks like maybe the `csne()` function would need more iterations to converge (?). I’ll have to read the reference they provide to better understand the method.