# Solve function for non full rank systems of equations

**URL:** https://discourse.julialang.org/t/solve-function-for-non-full-rank-systems-of-equations/54611
**Category:** Numerics
**Tags:** question, linearalgebra
**Created:** [February 4, 2021, 10:58am UTC](https://discourse.julialang.org/t/solve-function-for-non-full-rank-systems-of-equations/54611 "2021-02-04T10:58:10Z")
**Posts on this page:** 5
**Page:** 1

<div class="post-metadata">

### Author: ![Andres\_Legarra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andres_legarra/32/8786_2.png) [@Andres\_Legarra](https://discourse.julialang.org/u/Andres_Legarra)
#### Post date: [February 4, 2021, 10:58am UTC](https://discourse.julialang.org/t/solve-function-for-non-full-rank-systems-of-equations/54611/1 "2021-02-04T10:58:10Z")

</div>

Typical OLS ( X'X a = X'y ) with 2 or more dummy (or cross-classified) effects are not full rank and they have infinite solutions. (Minimal example below.) They are conceptually solved using generalized inverses, to obtain e.g. estimable functions with well known properties (e.g. Searle’s “Linear Models”). I am only interested in solving, so I don’t always want to use [pinv()](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.pinv) to explicitly get the generalized inverse. I don’t want to find singularities either - not useful and numerically unstable. Neither I want to explicitly use a package for linear models such as [GLM.jl](https://github.com/JuliaStats/GLM.jl) (Actually I want this for more than OLS.)

I had trouble finding appropriate solvers for non-full rank matrices in Julia; I tried `\` but it does not work because of singularity. Only `bicglstab` and `minres` from [IterativeSolvers.jl](https://julialinearalgebra.github.io/IterativeSolvers.jl/dev/linear_systems/minres/#IterativeSolvers.minres!) seem to work. Am I correct on using these for semi-positive definite X'X?

Also, am I missing something or is that all there is?

Thanks, Andres

```julia
#non full rank julia
# 2 effects perfectly balanced
X=[1 0 1 0; 0 1 0 1; 1 0 0 1; 0 1 1 0] .* 1.0
# 
y=[1 ; 2; 3 ; 4] .* 1.0

# sol= X'X \ X'y 
#does not work

using IterativeSolvers
sol = bicgstabl(X'X,X'y)
sol = minres(X'X,X'y)
# moore-penrose generalized inverse solution
sol = pinv(X'X)*X'y

```

---

<div class="post-metadata">

### Author: ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [February 4, 2021, 7:46pm UTC](https://discourse.julialang.org/t/solve-function-for-non-full-rank-systems-of-equations/54611/2 "2021-02-04T19:46:33Z")

</div>

How about `svd(X) \ y`?

---

<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: [February 4, 2021, 7:52pm UTC](https://discourse.julialang.org/t/solve-function-for-non-full-rank-systems-of-equations/54611/3 "2021-02-04T19:52:31Z")

</div>

You should generally use `qr(X, Val(true)) \ y` to get a least-square solution (defaulting to the minimum-norm solution if it is not unique, IIRC). For non-square `X`, it is equivalent to simply do `X \ y`.

You almost certainly do **not** want to form `X'X` **at all** if it is nearly singular (i.e. if `X` is badly conditioned), although this is deceptively the way least-squares problems are usually taught. See also [this post](https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/33).

---

<div class="post-metadata">

### Author: ![Andres\_Legarra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andres_legarra/32/8786_2.png) [@Andres\_Legarra](https://discourse.julialang.org/u/Andres_Legarra)
#### Post date: [February 5, 2021, 12:48pm UTC](https://discourse.julialang.org/t/solve-function-for-non-full-rank-systems-of-equations/54611/4 "2021-02-05T12:48:51Z")

</div>

Thanks. The OLS problem was a motivating example. I work with complex Mixed Models, never full rank, that we customarily solve using sparse matrices and iterative methods. Anyway IterativeSolvers seems to work fine to me.

---

<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: [February 5, 2021, 1:27pm UTC](https://discourse.julialang.org/t/solve-function-for-non-full-rank-systems-of-equations/54611/5 "2021-02-05T13:27:05Z")

</div>

> [@Andres\_Legarra](#):
>
> I work with complex Mixed Models, never full rank, that we customarily solve using sparse matrices and iterative methods.

There is also sparse QR. In general, I would look at sparse-direct methods first for sparse matrices, and only switch to iterative methods when the matrices get too large for sparse direct.)
