# Multiplication of submatrices problem

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

I think you have a couple surplus transposes, here

``````println(".  .  \$( L[k:n,1:r] * L[k,1:r]')")
``````

and here

``````L[k:n,r+1] = A[k:n,k] - ( L[k:n,1:r] * L[k,1:r]')
``````

Remove these two `'`s and it should work.

Seems to work, thank you.

Now I am really confused, however:
`A * B` should only work if the matrix sizes are consistent, e.g.,
A of size `M x K`, B of size `K x N` to yield a matrix of size `M x N`

`L[k:n,1:r]` should be size `(N-K+1) x R`, and
`L[k, 1:r]` should be size `1 x R`
How can they be multiplied with * without taking the transpose of the second factor?

I suspect that `L[k,1:r]` in fact is an `r` vector, i.e., 1-dimensional. And I think that in a Matrix times Vector product, the Vector is interpreted as a column matrix.

Darn, that would be inconsistent!
I could see this forcing if statementsâ€¦
I already have one in the code that should not be necessary

If you want to make matrices instead of vectors, you can index with a â€śtrivialâ€ť range, like `k:k` instead of `k`. See the difference in the types here:

``````julia> A = rand(4,4)
4Ă—4 Array{Float64,2}:
0.0250438   0.489103   0.631932  0.149773
0.389214    0.502628   0.882091  0.0734968
0.00247854  0.523021   0.892671  0.407138
0.895222    0.0263725  0.844777  0.206759

julia> A[2:2,3:4]
1Ă—2 Array{Float64,2}:
0.882091  0.0734968

julia> A[3:4,2:2]
2Ă—1 Array{Float64,2}:
0.5230213816602916
0.026372524205435566

julia> A[2,3:4]
2-element Array{Float64,1}:
0.8820909936504897
0.07349682319837192

julia> A[3:4,2]
2-element Array{Float64,1}:
0.5230213816602916
0.026372524205435566
``````

Great, thank you!

I do think this convention is a bit of a kludge, though!

The reasoning behind it is that if you want type stability, indexing a range needs to return the same type no matter what the range is. Thus range indexing in a dimension preserves that dimension. However scalar indexing of list should return an element. Thus a scalar index of an array should, by analogy drop that dimension.

3 Likes

So, given a matrix A,
`A[i,j]` is an entry in the matrix
`A[i,j:j]` is a vector (not a matrix), and
`A[i:i,j:j]` is a (sub) matrix?

Thank you for the explanation. I had missed it so far.

2 Likes