# Multiplication of submatrices problem

**URL:** <https://discourse.julialang.org/t/multiplication-of-submatrices-problem/33785>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [January 25, 2020, 7:24pm UTC](https://discourse.julialang.org/t/multiplication-of-submatrices-problem/33785 "2020-01-25T19:24:25Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![ea42gh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ea42gh/32/12000_2.png) [@ea42gh](https://discourse.julialang.org/u/ea42gh)\
**Post date:** [January 25, 2020, 7:24pm UTC](https://discourse.julialang.org/t/multiplication-of-submatrices-problem/33785/1 "2020-01-25T19:24:26Z")

</div>

I was trying to transcribe the algorithm in the appendix of [geninv](https://arxiv.org/pdf/0804.4809.pdf)

```julia
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

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

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [January 25, 2020, 8:58pm UTC](https://discourse.julialang.org/t/multiplication-of-submatrices-problem/33785/2 "2020-01-25T20:58:36Z")

</div>

I think you have a couple surplus transposes, here

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

```

and here

```julia
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.

---

<div class="post-metadata">

**Author:** ![ea42gh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ea42gh/32/12000_2.png) [@ea42gh](https://discourse.julialang.org/u/ea42gh)\
**Post date:** [January 25, 2020, 10:51pm UTC](https://discourse.julialang.org/t/multiplication-of-submatrices-problem/33785/3 "2020-01-25T22:51:11Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [January 25, 2020, 11:38pm UTC](https://discourse.julialang.org/t/multiplication-of-submatrices-problem/33785/4 "2020-01-25T23:38:49Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ea42gh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ea42gh/32/12000_2.png) [@ea42gh](https://discourse.julialang.org/u/ea42gh)\
**Post date:** [January 25, 2020, 11:42pm UTC](https://discourse.julialang.org/t/multiplication-of-submatrices-problem/33785/5 "2020-01-25T23:42:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [January 26, 2020, 1:04am UTC](https://discourse.julialang.org/t/multiplication-of-submatrices-problem/33785/6 "2020-01-26T01:04:56Z")

</div>

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

```

---

<div class="post-metadata">

**Author:** ![ea42gh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ea42gh/32/12000_2.png) [@ea42gh](https://discourse.julialang.org/u/ea42gh)\
**Post date:** [January 26, 2020, 1:13am UTC](https://discourse.julialang.org/t/multiplication-of-submatrices-problem/33785/7 "2020-01-26T01:13:24Z")

</div>

Great, thank you!

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 26, 2020, 1:56am UTC](https://discourse.julialang.org/t/multiplication-of-submatrices-problem/33785/8 "2020-01-26T01:56:20Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ea42gh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ea42gh/32/12000_2.png) [@ea42gh](https://discourse.julialang.org/u/ea42gh)\
**Post date:** [January 26, 2020, 2:03am UTC](https://discourse.julialang.org/t/multiplication-of-submatrices-problem/33785/9 "2020-01-26T02:03:22Z")

</div>

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.