Constructing a triangular power matrix from random variable


I would like to construct a square matrix which has the following form:

1.0, A, A^2, A^3, … A^n
0.0, 1.0, A, A^2, A^3, … A^n-1
0.0, 0.0, 1.0, A, A^2, … A^n-2

0.0, 0.0, 0.0, … 1.0

ie, the leading diagonal has ones, left/lower triangle is zeros, and right/upper triangle has power series, where ‘A’ is some random variable, eg pm.Normal()

the dimension of the matrix is determined by the input data.

how should I go about doing this ?


You can use an elementwise exponentiation to do this. See if this code helps:

import pymc3 as pm
import theano.tensor as tt
import numpy as np


# Dimension of square matrix
d = 5

C = np.zeros([d,d])

# Build a matrix in Numpy ahead of time
for i in range(d):
    C[i, i:] = np.arange(1, (d-i)+1)
with pm.Model() as model:
    A      = pm.Normal('A')
    matrix = pm.Deterministic('matrix',tt.triu(A ** pm.math.constant(C)))
    trace  = pm.sample()

yields the following output

array([[-5.25e-02,  2.76e-03, -1.45e-04,  7.61e-06, -4.00e-07],
       [ 0.00e+00, -5.25e-02,  2.76e-03, -1.45e-04,  7.61e-06],
       [ 0.00e+00,  0.00e+00, -5.25e-02,  2.76e-03, -1.45e-04],
       [ 0.00e+00,  0.00e+00,  0.00e+00, -5.25e-02,  2.76e-03],
       [ 0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00, -5.25e-02]])

thank you