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 :roll_eyes:

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