I was trying to transcribe the algorithm in the appendix of geninv
function pseudoinverse(G) # need to restrict types here
m,n = size(G)
if m < n
transpose=true
A=G*G'
n=m
else
transpose=false
A=G'G
end
dA = diag(A)
tol = 1e-9*minimum( dA[dA .> 0.] ) # catch dA=zeros()
L = fill( zero(A[1,1]), size(A) )
r = 0
for k = 1:n
if r > 0
println("indices: L[$k:$n,$r+1] = A[$k:$n,$k] - ( L[$k:$n,1:$r] * L[$k,1:$r]')")
println(". shape L[ $(size(L[k:n,1:r])) ]: coded as L[$k:$n,1:$r]")
println(". shape L[ $(size(L[k,1:r]')) ]: coded as L[$k,1:$r]'")
println(". . $( L[k:n,1:r] * L[k,1:r]')")
L[k:n,r+1] = A[k:n,k] - ( L[k:n,1:r] * L[k,1:r]')
else
L[k:n,r+1] = A[k:n,k]
end
r += 1
if L[k,r] > tol
L[k,r] = sqrt( L[k,r] )
if k < n
L[(k+1):n,r] = L[(k+1):n,r] ./ L[k,r]
end
else
r -= 1
end
end
L = L[:,1:r]
M = inv(L'L)
if transpose
G'*L*M*M*L'
else
L*M*M*L'*G'
end
end
pseudoinverse([1. 2 1; 2 2 1; 0 0 0])
It fails in the loop with
indices: L[3:3,2+1] = A[3:3,3] - ( L[3:3,1:2] * L[3,1:2]')
. shape L[ (1, 2) ]: coded as L[3:3,1:2]
. shape L[ (1, 2) ]: coded as L[3,1:2]'
which leaves me puzzled: the matrix sizes are consistent,
but somehow the transpose appears to be lost?!