Not able to multiply matrix with a vector of random Variable.

[prantik1998]

Hello everyone.

I am not able to multiply matrix with a vector of random variable

This is the code

function f(x,a,b,c,d)
    A = [a b; c d]
    A*x
end

@model function coin_toss_model(A,m,n)
    lambda = 0.5
    x = randomvar(m)
    y = datavar(Int64,n)
    x .~ Beta(1,lambda)
    u ~ f(x, A[1, 1], A[1, 2], A[2, 1], A[2, 2])
    for i in 1:m
        y[i] ~ Bernoulli(u[i])
    end
end

[HoangMHNguyen]

Hi @prantik1998! Thank you for playing around with RxInfer.jl. I think this issue isn't related to multiplying a matrix with a vector of random variables but it comes from your model definition. It seems that you've defined an invalid model. For example, Bernoulli only accepts a real number as input, not a vector.
Besides, the broadcasting feature isn't available in Rxinfer.jl yet.

[Nimrais]

@prantik1998 Hi, thank you for your question. Could you please provide more details about the probabilistic model you're trying to implement?

For instance, I'm not entirely clear on why it would be algebraically feasible to multiply your x by the matrix A. You seem to be attempting to multiply a 2x2 matrix with an m-dimensional vector. This operation can be not feasible due to the dimensional mismatch. Your model is valid only when m=2.

[Maximesncf]

Hi, I’m facing a similar issue in a model I’m trying to build.

It’s a kind of double HMM where the observations [0 or 1] depend on two hidden variables : one has two states and the other has three states. So I’m trying to do something like this:

function emission(z_1, z_2, B, t)
    emission_probs = [sum(z_1[t][i] * z_2[t][j] * B[i, j, k] for i in 1:2, j in 1:3) for k in 1:2]
    return emission_probs
end

Of course, this doesn’t work, and I don’t see how to perform the calculation...

@model function D_HMM(y,B)
    T_t ~ MatrixDirichlet(ones(2,2)) 
    T_b ~ MatrixDirichlet(ones(3,3)) 

    z_1_0 ~ Categorical([0.5, 0.5])  
    z_2_0 ~ Categorical([1/3, 1/3, 1/3])
    z_1_prev = z_1_0
    z_2_prev = z_2_0

    for t in eachindex(y)
        z_1[t] ~ Transition(z_1_prev, T_t)
        z_2[t] ~ Transition(z_2_prev, T_b)

        y := emission(z_1, z_2, B, t) # ???
 
        z_1_prev=z_1[t]
        z_2_prev=z_2[t]
    end
end

Any idea? Thank you in advance for your help.

[albertpod]

@Nimrais I think the only way is CVIProjection here?

P.S. @wouterwln isn't Transition node deprecated?

[Nimrais]

@albertpod I am not sure, kinda. The problem is that we would need to project on joint marginals over z_1 and z_2; there is no way right now in RxInfer to put such a projection constraint. So basically, because there is no edge z_1 in the model. Maybe I am missing smt.

[wouterwln]

I think you can get away with the general Transition node. As @albertpod noted, this will be deprecated in the newest version of RxInfer and renamed to DiscreteTransition, which will probably release somewhere next week. For now, you'll probably get away with the following:

@model function dual_hmm(y, B)
    T_t ~ MatrixDirichlet(ones(2,2)) 
    T_b ~ MatrixDirichlet(ones(3,3)) 

    z_1_0 ~ Categorical([0.5, 0.5])  
    z_2_0 ~ Categorical([1/3, 1/3, 1/3])
    z_1_prev = z_1_0
    z_2_prev = z_2_0

    for t in eachindex(y)
        z_1[t] ~ Transition(z_1_prev, T_t)
        z_2[t] ~ Transition(z_2_prev, T_b)

        y[t] ~ Transition(z_1[t], B, z_2[t])
 
        z_1_prev=z_1[t]
        z_2_prev=z_2[t]
    end
end

Or if you want to learn B:

@model function dual_hmm(y)
    B ~ TensorDirichlet(ones(2, 2, 3)) #(n_outputs x n_z_1 x n_z2)
    T_t ~ MatrixDirichlet(ones(2,2)) 
    T_b ~ MatrixDirichlet(ones(3,3)) 

    z_1_0 ~ Categorical([0.5, 0.5])  
    z_2_0 ~ Categorical([1/3, 1/3, 1/3])
    z_1_prev = z_1_0
    z_2_prev = z_2_0

    for t in eachindex(y)
        z_1[t] ~ Transition(z_1_prev, T_t)
        z_2[t] ~ Transition(z_2_prev, T_b)

        y[t] ~ Transition(z_1[t], B, z_2[t])
 
        z_1_prev=z_1[t]
        z_2_prev=z_2[t]
    end
end

For reference, in the newest RxInfer version, the code will change to:

@model function dual_hmm(y)
    B ~ DirichletCollection(ones(2, 2, 3)) #(n_outputs x n_z_1 x n_z2)
    T_t ~ DirichletCollection(ones(2,2)) 
    T_b ~ DirichletCollection(ones(3,3)) 

    z_1_0 ~ Categorical([0.5, 0.5])  
    z_2_0 ~ Categorical([1/3, 1/3, 1/3])
    z_1_prev = z_1_0
    z_2_prev = z_2_0

    for t in eachindex(y)
        z_1[t] ~ DiscreteTransition(z_1_prev, T_t)
        z_2[t] ~ DiscreteTransition(z_2_prev, T_b)

        y[t] ~ DiscreteTransition(z_1[t], B, z_2[t])
 
        z_1_prev=z_1[t]
        z_2_prev=z_2[t]
    end
end

Does this answer your question?

[wouterwln]

@Maximesncf Adding to my previous comment: It is wise to put at least some prior on B (or if you know it, disregard this comment) or on the transition matrices for z_1 and z_2. If you initialize your priors like I did (fully flat), your model has no way of knowing what caused the data and can hallucinate whatever latent representation for the hidden states. If you don't fully understand, I invite you to change the prior on the observation matrix in the HMM example on the RxInfer documentation (https://reactivebayes.github.io/RxInferExamples.jl/dev/examples/Basic%20Examples/Hidden%20Markov%20Model/Hidden%20Markov%20Model/) to a flat prior, and you will see that your model will start to hallucinate and the inference won't make sense. This is because without the slightly biased prior, your model has no incentive to trust any of its senses. I hope this additional comment helped you. Also, if you learn B, don't forget to add factorization constraints!

P.S. The rules for the Transition node are a bit messy, so I think you will have to wrap your data in the UnfactorizedData struct (as is shown here https://reactivebayes.github.io/RxInferExamples.jl/dev/examples/Basic%20Examples/Predicting%20Bike%20Rental%20Demand/Predicting%20Bike%20Rental%20Demand/) in order to make inference work. This is because the generalization of the Transition node only works when we learn q(y, z_1, z2)q(B), and when y is data, RxInfer will automatically try to learn q(y)q(z_1, z_2)1(B), which I haven't written any rules for in RxInfer (yet).

[Maximesncf]

Thank you for your response.

Indeed, I started by trying to work with the transition function (following your idea from here ) but without real success…

I get the error:
ERROR: Transition has 2 input interfaces (:in, :a), but 3 arguments provided.

T_t= [ 0.995 0.005;
          0     1]
T_b = [ 0.849 0.147 0.004;
          0.5   0.5   0;
          0     0     1]

B = [ 0.987 0.245 0.1; 
        0.1 0.75 0.9  ;;; 
        0 0 0          ; 
        1 1 1 ]

@model function dual_hmm(y, B, T_t, T_b)
    z_1_0 ~ Categorical([0.5, 0.5])  
    z_2_0 ~ Categorical([1/3, 1/3, 1/3])
    z_1_prev = z_1_0
    z_2_prev = z_2_0

    for t in eachindex(y)
        z_1[t] ~ Transition(z_1_prev, T_t)
        z_2[t] ~ Transition(z_2_prev, T_b)

        y[t] ~ Transition(z_1[t], B, z_2[t])
 
        z_1_prev=z_1[t]
        z_2_prev=z_2[t]
    end
end

constraints = @constraints begin
    q(z_1_0, z_1, z_2_0, z_2) = q(z_1_0, z_1)q(z_2_0,z_2)
end

returnvars = (
    z_1 = KeepLast(),
    z_2 = KeepLast(),
)

marginals = @initialization begin
    q(z_1) = vague(Categorical, 2)
    q(z_2) = vague(Categorical, 3)
end

result = infer(
    model = dual_hmm(B=B, T_b =T_b , T_t=T_t),
    data           = ( y = UnfactorizedData([ 1.0, 0.0, 0.0, 1.0 ]),),
    initialization = marginals,
    constraints    = constraints,
    returnvars     = returnvars,
    iterations     = 20,
    free_energy    = true
)

I mostly tried without wanting to learn the matrices, because I already know them to some extent… only the hidden states interest me here 😉.
That’s why I tried to manually create a sort of deterministic node, as I think the transition function isn’t quite suited to my use case (is it !?).

So, to try to build the node manually, I thought of retrieving the states to multiply them by the observation matrix. However, I started sinking into the abyss, trying several methods, with growing despair as my failure blossomed and my chances of success drifted further away—as Shakespeare might have put it :

y[t] := emission(z_1, z_2, B, t) 

or

 y[t] ~ emission(z_1, z_2, B, t)

ERROR: Delta node `emission` requires meta specification with the `where { meta = ... }` in the `@model` macro or with the separate `@meta` specification. See documentation for the `DeltaMeta`.

y[t] := Categorical([sum(z_1[t][i] * z_2[t][j] * B[i, j, k] for i in 1:2, j in 1:3) for k in 1:2])

ERROR: Indexing a single node label `z_1` with an index `[1]` is not allowed.

obs = emission(z_1, z_2, B, t)
y[t]  ~ Categorical(obs)

ERROR: MethodError: no method matching getindex(::GraphPPL.NodeLabel, ::Int64, ::Int64, ::Int64)
The function `getindex` exists, but no method is defined for this combination of argument types.

obs := emission(z_1, z_2, B, t)
y[t]  ~ Categorical(obs)

ERROR: MethodError: no method matching factornode(::ReactiveMP.UndefinedNodeFunctionalForm, ::Type{…}, ::Vector{…}, ::NTuple{…})
The function `factornode` exists, but no method is defined for this combination of argument types.`Type{TensorDirichlet}` has been used but the `ReactiveMP` backend does not support `Type{TensorDirichlet}` as a factor node.

And so on .... I'll spare you the rest, which I'm not very proud of.

Obviously, there’s something I’m missing

Do you think it would be more appropriate to try creating something with a @node instead? This seems feasible as long as it’s correct (or more in line with the philosophy of your library) to retrieve the hidden states as I’m trying to do (with the indexes)?

if I'm not mistaken, I'm trying to get this : $emission\_{probs_k}(t) = \sum_{i=1}^{2} \sum_{j=1}^{3} z_1[t][i] \cdot z_2[t][j] \cdot B[i, j, k] \quad \text{for} \quad k = 1, 2$

Finally, regarding your comment on initialization, I do have priors that could help. I didn’t think it was necessary to present all the code here. (Just to confirm, does it change anything about my "dual-HMM observation node issue" whether or not I need to learn T_b, T_t, and B?).

Thank you in advance for your help.

PS : I didn’t fully understand your remarks on the CVI projection, but should I understand that it’s not (yet) possible to manipulate/combine nodes this way?