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

```
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:

```
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.