How to solve for a transition matrix in Python?

Question

I have n vectors v_i, 1<=i<=n, and a starting vector v_0. I would like to solve for an nxn transition matrix A, s.t. A^i * v_0 = v_i for all i. What is the easiest way to implement this in python?

Answer 1

If v_0 ... v_n are linearly independent, then A = [v_1 .. v_n] * [v_0 ... v_{n-1}]^-1, i.e. A = np.hstack((v_1, ..., v_n)) @ inv(np.hstack((v_0, ..., v_n_1))).

Otherwise, no unique solution.

Answer 2

This is a pretty tough question since you have a nonlinear set of equations. I've attempted to solve it numerically using scipy.optimize.fsolve, but I can only sometimes get reliable solutions once I'm up to n = 3.

import scipy.optimize as sco
import numpy as np

# Set up problem
n = 3

# Generate a random transition array and v_0
real_A = np.random.randint(0, 10, (n, n))
v_i = {'v0' : np.random.randint(0, 10, n)}

# Calculate v_1 to v_n, all v_i stored in dictionary
for i in range(1, n+1):
    v_i[f'v{i}'] = real_A @ v_i[f'v{i-1}']

# Create function where retrieve_A(real_A, v_i) returns all zeros
def retrieve_A(x, v_i):
    n = len(v_i) - 1
    assert len(x) == n*n, 'Not enough arguments in x'
    
    A = np.reshape(x, (n, n))
    
    
    out = {'v0': v_i['v0']}
    
    # These are "v_i" generated by our transition matrix
    for i in range(1, n+1):
        out[f'v{i}'] = A @ out[f'v{i-1}']
    
    # Return comparison of "out_i" to our actual "v_i"
    return np.array([out[f'v{i}'] - v_i[f'v{i}'] for i in range (1, n+1)]).ravel()


# Use scipy.optimize.fsolve to solve the equation set by retrieve_A and hopefully get our A matrix
retrieved_A = sco.fsolve(retrieve_A, [1]*(n*n), args = (v_i)).reshape(n, n)


print('Real A:', real_A)
print('Solved A:', retrieved_A)