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