Solve function for non full rank systems of equations

Specific Domains Numerics
,
1

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() 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 (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 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

#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
2

How about svd(X) \ y?

3

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.

3 Likes
4

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.

5

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

1 Like