Calculation of quadratic form using broadcasting in Julia

Question

I want to calculate a vector of a quadratic form, extracting the submatrix from 3 by 3 by 5 arrays. However, I cannot make the quadratic form using broadcasting (i.e., macro "@."). When using “for” statement, we can calculate the vector of the quadratic form. I have no idea how to conduct matrix operations using “@.” (I am reluctant to expand the quadratic form to calculate the vector.)

By contrast, the inner product is computable using “@.”.

The example code is as follows:

using LinearAlgebra

a1=[5 7 2; 2 1 5; 6 2 3]
a2=[2 7 1; 3 7 2; 1 2 3]
a3=[8 5 9; 1 1 3; 2 2 3]
a4=[2 5 6; 3 5 1; 1 1 1]
a5=[7 8 1; 5 1 3; 1 5 2]
z=cat(a1,a2,a3,a4,a5,dims=3)


##### case of inner product
x=zeros(5,3)
wz = reshape([],0)
for k in 1:5
    w = hcat(z[[1],[1],k], z[2,2,k]) * hcat(z[[1],[1],k], z[[2],[2],k])'
    #println(w)
    wz=vcat(wz, w)
end
@. wz=convert(Float64,wz)
wz=Matrix{Float64}(wz)
x[:,3]=wz
# [inner product] same result, the 3rd column vector [26.0, 53.0, 65.0, 29.0, 50.0]
display(x)

x=zeros(5,3)
@. x[:,3] = dot(hcat(z[1,1,:],z[2,2,:]), hcat(z[1,1,:],z[2,2,:])) # ok, working
# [inner product] same result, the 3rd column vector [26.0, 53.0, 65.0, 29.0, 50.0]
display(x) 


##### case of quadratic form
x=zeros(5,3)
wy = reshape([],0)
for k in 1:5
    w = hcat(z[[1],[1],k], z[[2],[2],k]) * z[[1,3],[1,3],k] * hcat(z[[1],[1],k], z[[2],[2],k])'
    #println(w)
    wy=vcat(wy, w)
end
@. wy=convert(Float64,wy)
wy=Matrix{Float64}(wy)
x[:,3]=wy
# [quadratic form] distinct result, the 3rd column vector [168.0, 183.0, 603.0, 103.0, 359.0]
display(x)

# generating five 2 by 2 matrices, distinct result
@. dot(hcat(z[[1],[1],:],z[[2],[2],:]), z[[1,3],[1,3],:], hcat(z[[1],[1],:],z[[2],[2],:]))

# obtaining ERROR: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 2 and 5")
@. dot(hcat(z[1,1,:],z[2,2,:]), z[[1,3],[1,3],:], hcat(z[1,1,:],z[2,2,:]))

Would you mind giving helps and suggestions how to get the calculation of 3rd column vector [168.0, 183.0, 603.0, 103.0, 359.0] (which is made from the quadratic form) in the above code using "@."?

Answer 1

EDIT:

Perhaps the question is about specifically how to make broadcasting work in this case. If so:

@views dot.(vcat.(z[1,1,:],z[2,2,:]),getindex.(Ref(z),Ref([1,3]),Ref([1,3]),axes(z,3)),vcat.(z[1,1,:],z[2,2,:]))

should be a possible clarification. Or with the @. macro (though it doesn't seem simpler):

@. dot(vcat(z[1,1,:],z[2,2,:]),getindex($Ref(z),$Ref([1,3]),$Ref([1,3]),$axes(z,3)),vcat(z[1,1,:],z[2,2,:]))

ORIGINAL:

One way to calculate this:

[
  [z[1,1,k] z[2,2,k]]*z[[1,3],[1,3],k]*[z[1,1,k] z[2,2,k]]' |> first 
  for k ∈ axes(z,3)
]

giving:

5-element Vector{Int64}:
 168
 183
 603
 103
 359

(the |> first turns 1x1 matrix into scalar)

Option 2:

[let t = z[[1,3],[1,3],k] ; sum(z[i,i,k]*t[i,j]*z[j,j,k] for i ∈ (1,2), j ∈ (1,2)) ; end for k ∈ 1:5]

or:

[let t = z[[1,3],[1,3],k], v = [z[1,1,k],z[2,2,k]] ; dot(v,t,v) ; end for k ∈ 1:5]

or (this is pretty cool):

map((z;t=z[[1,3],[1,3]],v=[z[1,1],z[2,2]])->dot(v,t,v), eachslice(z,dims=3))