Outer product matrix multiplication

user_19244821 0
Viewed 99

Trying to find the outer product of two matrices A and B

Here is what I have attempted:

function product(A, B)
    n_a, m_a = size(A)
    n_b, m_b = size(B)
    AB = Array{Float64}(undef, n_a, m_b)
    for i in 1:n_a
        for j in 1:m_b
            AB[i, j] = A[:, j] * B[j, :]'
        end
    end
    return AB
end

product(A, B)

I get an error when attempting to run: Cannot `convert` an object of type Matrix{Float64} to an object of type Float64

julia matrix-multiplication
3 Answers

I'm not exactly sure what you mean by "outer product of matrices", I'm only familiar with an outer product of vectors (which creates a matrix). Could you clarify what output you're looking for with an example?

In any event to address the immediate issue: The stacktrace is pointing to this line in your function:

AB[i, j] = A[:, j] * B[j, :]'

Let's take two example matrices and see what happens for i=j=1:

julia> x = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> y = [2 3; 4 5]
2×2 Matrix{Int64}:
 2  3
 4  5

julia> x[:, 1] * y[1, :]'
2×2 Matrix{Int64}:
 2  3
 6  9

so the way you are slicing and transposing your matrices means you are calculating an outer product (as I know it) of two vectors, which gives you a matrix. Given that AB is defined as a matrix of Float64s, the location AB[i, j] can only hold a Float64 value, but you are trying to assign a Matrix{Float64} to it.

Again I'm not sure what exactly you are trying to achieve here, but if it's just a "normal" matrix multiplication, you should have

AB[i, j] = A[i, :]' * B[:, j]

in your inner loop. Making that change gives me:

julia> product(x, y) == x * y
true

for the x and y defined above.

user_24998920

Based on the image you linked to in the comment, you're calculating an array of matrices. This can be done as follows:

AB = [A[:,j]*B[j,:]' for j=1:size(A,2)]

Though, I don't believe that's the correct definition for matrix outer product. I think it should be like this:

AB = [r*c' for r in eachrow(A), c in eachcol(B)]

This will give you a [n_a*m_b] matrix of [m_a*n_b] matrices.

user_39431700

From the linked image and the text, I think the formula tries to show another representation for the usual matrix product, which is useful, especially when doing rank decomposition of matrices (as sums of rank one matrices). As a concrete example, first defining a couple of matrices:

julia> using Random
julia> Random.seed!(1234);

julia> A = rand(1:10,(2,3))
2×3 Matrix{Int64}:
 4  3  4
 6  9  4

julia> B = rand(1:10,(3,4))
3×4 Matrix{Int64}:
 10  8  1   7
  8  6  2  10
  5  8  5   7

Now define a product as a sum of outer products of vectors, and see it gives the same result as a usual matrix product:

julia> product(A, B) = sum( [ A[:,i] * B[i,:]' for i=1:size(A,2) ] )
product (generic function with 1 method)

julia> product(A,B)
2×4 Matrix{Int64}:
  84   82  30   86
 152  134  44  160

julia> A*B
2×4 Matrix{Int64}:
  84   82  30   86
 152  134  44  160

In this example, the resulting matrix can be at most rank 3, the number of columns in A (and rows in B). In fact, it is rank 2 because of number of rows, but in general this product is done on tall times flat matrices and then the rank restriction is meaningful.

user_35808700
Related