CODEWITHSUNDEEP
pythonlinear-algebramatrix-multiplicationpolynomial-math
Soumyadeep sarma
Soumyadeep sarma

Reputation: 11

Matrix multiplication on Python to find a polynomial

I wish to find a method to multiply matrices (with other matrices and vectors) whose elements contain polynomial expressions and to obtain the final results in polynomial format.

I have a 10 by 10 matrix A, whose elements are function of a variable z. Now, I wish to exponentiate this matrix to some power (say L), and then multiply it with a 10 by 1 column vector. What I need, is the polynomial expression of first element of the resulting final 10 by 1 column vector (say f(z)), as a function of z, and then find its roots (it will be a polynomial with degree L, so L roots). How can I implement this in Python? I have been doing it in Mathematica up until recently, where it is easy to multiply polynomial expressions, but I need to scale up the matrix sizes, and it is difficult to do so in Mathematica alone. Any such matrix in Python asks us to define the value of z first, so it will effectively only compute the value of f(z) at the initialized value of z, hence the issue.

Upvotes: 1

Views: 128

Answers (1)

tbrugere
tbrugere

Reputation: 1623

Sympy might be able to do that out of the box.

Otherwise with just numpy and scipy, here is a simple way to do this:

  1. Encode your matrix A as a matrix of padded vector 10 × 10 × max_degree * (L + 1) where eg. A[3, 2, 0] is the constant term factor of the polynomial at coordinates 3, 2.

Note: they should be padded to max_degree * (L + 1) for there to be enough coordinates even after the multiplications!

  1. apply fourier transform
A_fft = np.fft(A, axis=-1)
  1. Now remark that since polynomial multiplication is equivalent to Hadamard (value-wise) multiplication in the Fourier space, you can just use the following:
A_fft_power_i = A_fft

for i in range(1, L):
   # matrix multiplication of a matrix of vectors
   # where component-wise multiplication is used to multiply two vectors
   A_fft_power_i = np.einsum("ijc,jkc->ikc", A_fft_power_i, A)

A_fft_power_L = A_fft_power_i

#assuming v is the vector you want to multiply by 
#(and that it is a vector of polynomials too encoded the same as A 
# - if that isn't the case you can always consider that it is a vector of constant polynomials)
v_fft = np.fft(v, axis=-1)
A_fft_power_L_times_v = np.einsum("ijc,jc->ic", A_fft_power_L, v)

(note: this isn't very optimized... if computing time is of the essence, you could reimplement a version of fast exponentiation)

  1. Get back the non-fft value
result = np.ifft(A_fft_power_L_times_v[0])

Upvotes: 1

Related Questions