Skip to content

Documentation under construction

📝 Creating polynomials

In PolyAny a polynomial can be created in many forms. The most flexible way is defining a matrix of exponents and a vector of coefficients.

1⃣ From exponents and coefficients

The exponents matrix contains the exponents of each monomial of the polynomial, and the coefficients vector contains the corresponding scalars multipliers.

Basic example

Consider the monomial \(M(\mathbf{x}) = 5\,x_1^2\,x_2\,x_3^4\,x_5\), its exponents and coefficient are given by:

exponents = [[2, 1, 4, 0, 1]]  # (1)!
coefficient = [5]
  1. When declaring the exponents matrix, variables are ordered increasingly, i.e., \([x_1,\,x_2,\,x_3,\,\dots,\,x_n]\).

Note that the exponent matrix represents the complete monomial, meaning that it must include all variables, even those raised to the power of \(0\).

Example

If we want to represent a polynomial, we declare the exponents and coefficients of each monomial. Consider:

\[ P(\mathbf{x}) = 3\,x_1\,x_2 + 5\,x_1^2\,x_2\,x_3^4\,x_5 + 4\,x_4^4\,x_5^3 \]

The monomials of \(P(\mathbf{x})\) are:

  • \(M_1(\mathbf{x}) = 3\,x_1\,x_2\)
  • \(M_2(\mathbf{x}) = 4\,x_4^4\,x_5^3\)
  • \(M_3(\mathbf{x}) = 5\,x_1^2\,x_2\,x_3^4\,x_5\)

To represent \(P(\mathbf{x})\), we define:

exponents = [[1, 1, 0, 0, 0], [0, 0, 0, 4, 3], [2, 1, 4, 0, 1]]
coefficients = [3, 4, 5]

It is very important that the coefficients appear in the same order as their corresponding rows in the exponents matrix.

To create a polynomial from an exponents matrix and a coefficients vector in PolyAny, use the following syntax:

from polyany import Polynomial

exponents = [[0, 0], [1, 0], [0, 1]]  # (1)!
coefficients = [1, 2, 3]  # (2)!

poly = Polynomial(exponents, coefficients)
  1. The exponents matrix can also be a nested tuples or a NumPy 2D array
  2. The coefficients vector can also be a tuple or a NumPy 1D array

The code above creates the following multivariate polynomial:

\[ P(\mathbf{x}) = 1 + 2\,x_1 + 3\,x_2 \]

2⃣ Univariate Polynomials

Univariate Polynomials

A univariate polynomial \(P(x_1)\) is a polynomial that depends on a single variable \(x_1\). An example of univariate polynomial is:

\[ P(x_1) = 4 + 3 x_1 + 2 x_1^2 \]

To create univariate polynomials, a simpler syntax can be used, it only requires a coefficients vector.

from polyany import Polynomial

coefficients = [1, 5, 8, 9]
univar_poly = Polynomial.univariate(coefficients)

The coefficients are related with power of the variable \(x_1\) in increasing degree. The code above creates the polynomial:

\[ P(x_1) = 1 + 5\,x_1 + 8\,x_1^2 + 9\,x_1^3 \]

3⃣ Quadratic forms

Quadratic forms

A quadratic form is a second-degree homogeneous multivariate polynomial. It can be represented using a symmetric matrix \(A \in \mathbb{R}^{n \times n}\) as:

\[ P(\mathbf{x}) = \sum_{i=1}^{n}\sum_{j=1}^{n} a_{ij} x_i x_j = \mathbf{x}^{\top}A\mathbf{x} \]

For further reading, see Linear Algebra Done Right by Sheldon Axler

To create a quadratic form a square matrix must be provided.

from polyany import Polynomial

matrix = [[1, 2], [2, 3]]  # (1)!
quadratic_form_poly = Polynomial.quadratic_form(matrix)
  1. The input matrix can also be nested tuples or a NumPy 2D array

The code above creates the polynomial:

\[ P(\mathbf{x}) = x_1^2 + 4\,x_1\,x_2 + 3\,x_2^2 \]

Non-symmetric matrices

If the input matrix \(A\) is not symmetric, a warning is raised and it's symmetric part \(A_{\mathrm{sym}}\) is used instead, where

\[ A_{\mathrm{sym}} = \frac{1}{2} \left( A + A^{\top} \right) \]
Example

Suppose the input matrix is:

\[ A = \begin{bmatrix} 1 & 6 \\ 0 & 2 \end{bmatrix} \]

Its symmetric part is:

\[ A_{\mathrm{sym}} = \begin{bmatrix} 1 & 3 \\ 3 & 2 \end{bmatrix} \]

which corresponds to the polynomial:

\[ P(\mathbf{x}) = x_1^2 + 6\,x_1\,x_2 + 2\,x_2^2 \]

🔢 Evaluating polynomials

Polynomials can be easily evaluated in PolyAny by treating the polynomial object as a callable function. The input argument (point) must be a vector with n_vars components.

>>> poly = Polynomial.univariate([1, 2, 3])
>>> poly([2])  # (1)!
np.float64(17.0)
  1. Even for univariate polynomials, point needs to be a list, tuple or NumPy 1D-array.

For multivariate polynomials:

>>> matrix = [[1, 2, 3], [2, 4, 5], [3, 5, 6]]
>>> poly = Polynomial.quadratic_form(matrix)
>>> poly([0, 0, 0])
np.float64(0.0)
>>> poly([1, 1, 1])
np.float64(31.0)

🟰 Comparing polynomials

In PolyAny, Polynomial objects support equality comparisons (==) with other polynomials, but do not support ordering comparisons (<, <=, >, >=), which raises a TypeError.

Two polynomials are considered equal if and only if they have:

  • The same number of variables (n_vars attribute)
  • The same total degree (degree attribute)
  • The same coefficients (coefficients attribute)1

Comparison with other types (sequences, scalars, NumPy arrays) always returns NotImplemented.

Example

Internally, the polynomial object stores the coefficients and exponents in an ordered way, meaning that a polynomial is created regardless of the order of the input coefficients and exponents.

>>> poly1 = Polynomial([[0, 0], [1, 0], [0, 1]], [1, 2, 3])
>>> poly2 = Polynomial([[0, 1], [0, 0], [1, 0]], [3, 1, 2])
>>> poly1 == poly2
True

✂ Pruning

Pruning is the process of removing empty monomials of a polynomial. The Polynomial object stores the exponents and coefficients provided by the user in a ordered way.

>>> poly = Polynomial([[0, 0], [0, 1], [1, 0], [1, 1]], [1, 0, 2, 0])
>>> poly.exponents
array([[0, 0],
       [1, 0],
       [0, 1],
       [1, 1]])
>>> poly.coefficients
array([1., 2., 0., 0.])

When a polynomial is pruned, all empty monomials are removed, that is, the entries in exponents whose associated coefficients are exactly zero, which have no effect on the polynomial behavior.

To prune a polynomial, use the prune method:

>>> pruned = poly.prune()
>>> pruned.exponents
array([[0, 0],
       [1, 0]])
>>> pruned.coefficients
array([1., 2.])

The pruned polynomial retains only the first and second term, which are the non-empty monomials.

🗜 Squeezing

Squeezing is the process of removing extra variables of a polynomial. That is, it removes the columns with all zeros in exponents.

Note

The squeezing process doesn't alter coefficients and degree attributes.

>>> poly = Polynomial([[0, 0, 0], [0, 1, 0], [0, 2, 0]], [1, 2, 3])
>>> poly.exponents
array([[0, 0, 0],
       [0, 1, 0],
       [0, 2, 0]])

In the example above, the variables \(x_1\) and \(x_3\) are not utilized, so they can be removed by using the squeeze method:

>>> squeezed = poly.squeeze()
>>> squeezed.exponents
array([[0],
       [1],
       [2]])

➕ Addition and subtraction

In PolyAny, Polynomial objects can be added to or subtracted from scalars2 and other scalar polynomials.

>>> poly = Polynomial.univariate([1, -2, 3])
>>> poly
1 - 2*x_1 + 3*x_1^2
>>> poly + 5
6 - 2*x_1 + 3*x_1^2
>>> poly - 1
-2*x_1 + 3*x_1^2

For addition/subtraction between polynomials:

>>> another_poly = Polynomial([[0, 0], [1, 0], [0, 1], [1, 1]], [1, -2, 3, -4])
>>> another_poly
1 - 2*x_1 + 3*x_2 - 4*x_1*x_2
>>> poly + another_poly
2 - 4*x_1 + 3*x_1^2 + 3*x_2 - 4*x_1*x_2
>>> poly - another_poly
3*x_1^2 - 3*x_2 + 4*x_1*x_2

Similarly, MatrixPolynomial objects support addition and subtraction with scalars2, matrices3, and other matrix polynomials.

Interaction between scalar and matrix polynomials

MatrixPolynomial objects cannot operate with Polynomial objects.

>>> C_1 = np.eye(3)
>>> C_2 = np.ones((3, 3))
>>> C_3 = np.arange(9).reshape(3, 3)
>>> mpoly = MatrixPolynomial([[0], [1], [2]], [C_1, C_2, C_3])
>>> mpoly
[[1. 0. 0.]    [[1. 1. 1.]        [[0. 1. 2.]
 [0. 1. 0.]     [1. 1. 1.]         [3. 4. 5.]
 [0. 0. 1.]] +  [1. 1. 1.]]*x_1 +  [6. 7. 8.]]*x_1^2
>>> mpoly + 10  # (1)!
[[11. 10. 10.]    [[1. 1. 1.]        [[0. 1. 2.]
 [10. 11. 10.]     [1. 1. 1.]         [3. 4. 5.]
 [10. 10. 11.]] +  [1. 1. 1.]]*x_1 +  [6. 7. 8.]]*x_1^2
>>> mpoly + np.eye(3)  # (2)!
[[2. 0. 0.]    [[1. 1. 1.]        [[0. 1. 2.]
 [0. 2. 0.]     [1. 1. 1.]         [3. 4. 5.]
 [0. 0. 2.]] +  [1. 1. 1.]]*x_1 +  [6. 7. 8.]]*x_1^2
  1. Similarly to NumPy, operations are performed using broadcasting.
  2. You could use nested lists or nested tuples.

Operating between matrix polynomials:

>>> another_mpoly = MatrixPolynomial([[1]], np.ones((1, 3, 3)))
>>> mpoly - another_mpoly
[[1. 0. 0.]    [[0. 1. 2.]
 [0. 1. 0.]     [3. 4. 5.]
 [0. 0. 1.]] +  [6. 7. 8.]]*x_1^2

✖ Multiplication and division

Interaction between scalar polynomials and matrix polynomials

MatrixPolynomial objects cannot operate with Polynomial objects.

In PolyAny, Polynomial objects can be multiplied with other polynomials and scalars2.

>>> poly = Polynomial.univariate([10, -20, 5])
>>> poly
10 - 20*x_1 + 5*x_1^2
>>> poly * 2
20 - 40*x_1 + 10*x_1^2

Multiplying two polynomials:

>>> poly1 = Polynomial.univariate([1, -2, 3])
>>> poly1
1 - 2*x_1 + 3*x_1^2
>>> poly2 = Polynomial([[0, 0], [1, 1]], [3, 3])
>>> poly2
3 + 3*x_1*x_2
>>> poly1 * poly2
3 - 6*x_1 + 9*x_1^2 + 3*x_1*x_2 - 6*x_1^2*x_2 + 9*x_1^3*x_2

MatrixPolynomial objects support element-wise multiplication with scalars2, matrices3, and other matrix polynomials.

>>> mpoly = MatrixPolynomial([[0], [1]], [np.tri(3), np.vander([1, 2, 3])])
>>> mpoly
[[1. 0. 0.]    [[1. 1. 1.]
 [1. 1. 0.]     [4. 2. 1.]
 [1. 1. 1.]] +  [9. 3. 1.]]*x_1
>>> mpoly * 2 # (1)!
[[2. 0. 0.]    [[ 2.  2.  2.]
 [2. 2. 0.]     [ 8.  4.  2.]
 [2. 2. 2.]] +  [18.  6.  2.]]*x_1
>>> mpoly * np.eye(3) # (2)!
[[1. 0. 0.]    [[1. 0. 0.]
 [0. 1. 0.]     [0. 2. 0.]
 [0. 0. 1.]] +  [0. 0. 1.]]*x_1
  1. Similarly to NumPy, operations are performed using broadcasting.
  2. You could use nested lists or nested tuples.

Division between polynomials

Currently, division can only be performed between Polynomial objects and scalars2. In the future, it is possible that division between polynomial objects will be supported.

Dividing a Polynomial object by a scalar2:

>>> poly = Polynomial.univariate([10, -20, 5])
>>> poly
10 - 20*x_1 + 5*x_1^2
>>> poly / 5
2 - 4*x_1 + x_1^2

Dividing a MatrixPolynomial object by a scalar2:

>>> mpoly = MatrixPolynomial(
  [[0, 0], [1, 0], [0, 1]],
  [np.tri(3), np.eye(3), np.vander([3, 1, 4])],
)
>>> mpoly
[[1. 0. 0.]    [[1. 0. 0.]        [[ 9.  3.  1.]
 [1. 1. 0.]     [0. 1. 0.]         [ 1.  1.  1.]
 [1. 1. 1.]] +  [0. 0. 1.]]*x_1 +  [16.  4.  1.]]*x_2
>>> mpoly / 4
[[0.25 0.   0.  ]    [[0.25 0.   0.  ]        [[2.25 0.75 0.25]
 [0.25 0.25 0.  ]     [0.   0.25 0.  ]         [0.25 0.25 0.25]
 [0.25 0.25 0.25]] +  [0.   0.   0.25]]*x_1 +  [4.   1.   0.25]]*x_2

🌀 Matrix multiplication

In PolyAny, matrix multiplication can be performerd on MatrixPolynomial objects and matrices3.

Warning

Scalars2 are not accepted in matrix multiplication, use element-wise multiplication instead.

>>> mpoly = MatrixPolynomial([[0], [1]], [[[3, 1],[4, 1]], np.tri(2)])
>>> mpoly
[[3. 1.]    [[1. 0.]
 [4. 1.]] +  [1. 1.]]*x_1
>>> mpoly @ np.ones((2, 2))
[[4. 4.]    [[1. 1.]
 [5. 5.]] +  [2. 2.]]*x_1
>>> np.ones((2,2)) @ mpoly
[[7. 2.]    [[2. 1.]
 [7. 2.]] +  [2. 1.]]*x_1

Matrix product of two matrix polynomials:

>>> another_mpoly = MatrixPolynomial([[1, 0], [0, 1]], [np.vander([1, 2]), np.diag([1, 2])])
>>> another_mpoly
[[1. 1.]        [[1. 0.]
 [2. 1.]]*x_1 +  [0. 2.]]*x_2
>>> mpoly @ another_mpoly
[[5. 4.]        [[1. 1.]          [[3. 2.]        [[1. 0.]
 [6. 5.]]*x_1 +  [3. 2.]]*x_1^2 +  [4. 2.]]*x_2 +  [1. 2.]]*x_1*x_2
>>> another_mpoly @ mpoly
[[ 7.  2.]        [[2. 1.]          [[3. 1.]        [[1. 0.]
 [10.  3.]]*x_1 +  [3. 1.]]*x_1^2 +  [8. 2.]]*x_2 +  [2. 2.]]*x_1*x_2

Matrix multiplication is not-commutative

In the examples above, notice that operand_1 @ operand_2 is different of operand_2 @ operand_1.

➰ Partial derivatives

The partial derivatives of polynomials can be evaluated by using the partial method. Let's consider the polynomial:

\[ P(\mathbf{x}) = 10 + 2\,x_1^2\,x_2\,x_3 + 5\,x_1\,x_2^3\,x_3^2 \]

This can be declared in PolyAny as:

>>> poly = Polynomial([[0, 0, 0], [2, 1, 1], [1, 3, 2]], [10, 2, 5])

Their first partial derivatives are:

\[ \begin{cases} \displaystyle\frac{\partial P(\mathbf{x})}{\partial x_1} = 4\,x_1\,x_2\,x_3 + 5\,x_2^3\,x_3^2 \\[.5em] \displaystyle\frac{\partial P(\mathbf{x})}{\partial x_2} = 2\,x_1^2\,x_3 + 15\,x_1\,x_2^2\,x_3^2 \\[.5em] \displaystyle\frac{\partial P(\mathbf{x})}{\partial x_3} = 2\,x_1^2\,x_2 + 10\,x_1\,x_2^3\,x_3 \end{cases} \]

which can be obtained in PolyAny as:

>>> poly.partial(0)  # (1)!
4*x_1*x_2*x_3 + 5*x_2^3*x_3^2
>>> poly.partial(1)
2*x_1^2*x_3 + 15*x_1*x_2^2*x_3^2
>>> poly.partial(2)
2*x_1^2*x_2 + 10*x_1*x_2^3*x_3
  1. The method partial uses a zero-based index.

Matrix Polynomials

Definition

A matrix polynomial is a polynomial whose coefficients are matrices of the same shape. An example of a matrix polynomial is:

\[ P(\mathbf{x}) = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}\,x_1 + \begin{bmatrix} 3 & 1 \\ 4 & 5 \end{bmatrix}\,x_2 \]

In PolyAny, a matrix polynomial can be created from exponents and coefficients, similar to scalar polynomials, using the MatrixPolynomial class.

The shape attribute

The shape (number of rows and the number of columns) of the matrices in the polynomial can be obtained by the shape attribute. Following the example:

>>> mpoly.shape
(2, 2)

1⃣ From exponents and coefficients

If we want to declare the matrix polynomial in the definition above, we define:

from polyany import MatrixPolynomial

exponents = [[1, 0], [0, 1]]  # (1)!
C_1 = [[1, 0], [0, 1]]  # (2)!
C_2 = [[3, 1], [4, 5]]
coefficients = [C_1, C_2]  # (3)!

mpoly = MatrixPolynomial(exponents, coefficients)
  1. The exponents matrix can also be nested tuples or a NumPy 2D array.
  2. You could use np.eye method.
  3. The coefficients array can also be nested tuples or a NumPy 3D array.

2⃣ From a scalar polynomial

A Polynomial object can be converted into a MatrixPolynomial with the from_scalar classmethod.

The method recieves three attributes: the scalar polynomial to be converted, the desired shape of the resultant polynomial, and the conversion method.

Conversion methods

The default conversion method is ones, which uses a ones matrix to expand the scalar coefficients. The eye method uses an identity matrix instead.

>>> poly = Polynomial([[0, 0], [1, 0], [0, 1]], [1, 2, 3])
>>> poly
1 + 2*x_1 + 3*x_2
>>> MatrixPolynomial.from_scalar(poly, shape=(3, 2)) # (1)!
[[1. 1.]    [[2. 2.]        [[3. 3.]
 [1. 1.]     [2. 2.]         [3. 3.]
 [1. 1.]] +  [2. 2.]]*x_1 +  [3. 3.]]*x_2
>>> MatrixPolynomial.from_scalar(poly, shape=(2, 2), method="eye")
[[1. 0.]    [[2. 0.]        [[3. 0.]
 [0. 1.]] +  [0. 2.]]*x_1 +  [0. 3.]]*x_2
  1. By default, the conversion method is ones.

🇹 Transposition

To transpose a MatrixPolynomial object you can use the T property.

Example

Create a matrix polynomial:

>>> mpoly = MatrixPolynomial([[1, 0], [0, 1]], [np.tri(3), np.arange(9).reshape(3, 3)])
>>> mpoly
[[1. 0. 0.]        [[0. 1. 2.]
 [1. 1. 0.]         [3. 4. 5.]
 [1. 1. 1.]]*x_1 +  [6. 7. 8.]]*x_2

obtains its transpose with the T property:

>>> mpoly.T
[[1. 1. 1.]        [[0. 3. 6.]
 [0. 1. 1.]         [1. 4. 7.]
 [0. 0. 1.]]*x_1 +  [2. 5. 8.]]*x_2

🧩 Concatenation

Simple concatenation

Definition

Simple concatenation is the process of concatenating the coefficient matrices of several polynomials (vertically or horizontally) with respect to each monomial. If a monomial exists in one polynomial but not in the others, a zeros matrix of appropriate shape is utilized.

Vertical concatenation

Consider these matrix polynomials:

\[ P_1(\mathbf{x}) = \underbrace{ \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}}_{A_1}\,x_1\,x_2 + \underbrace{ \begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix}}_{B_1}\,x_2^2 \]
\[ P_2(\mathbf{x}) = \underbrace{ \begin{bmatrix} 0 & 1 \\ 2 & 3 \\ 4 & 5 \end{bmatrix}}_{A_2} + \underbrace{ \begin{bmatrix} 3 & 3 \\ 3 & 3 \\ 3 & 3 \end{bmatrix}}_{B_2}\,x_2^2 \]

Their vertical concatenation is:

\[ \begin{bmatrix} P_1(\mathbf{x}) \\ P_2(\mathbf{x}) \end{bmatrix} = \begin{bmatrix} \mathbf{0}_{2 \times 2} \\ A_2 \end{bmatrix} + \begin{bmatrix} A_1 \\ \mathbf{0}_{3 \times 2} \end{bmatrix}\,x_1\,x_2 + \begin{bmatrix} B_1 \\ B_2 \end{bmatrix}\,x_2^2 \]

which is:

\[ \begin{bmatrix} P_1(\mathbf{x}) \\ P_2(\mathbf{x}) \end{bmatrix} = \begin{bmatrix} 0 & 0 \\ 0 & 0 \\ \hline 0 & 1 \\ 2 & 3 \\ 4 & 5 \end{bmatrix} + \begin{bmatrix} 1 & 1 \\ 1 & 1 \\ \hline 0 & 0 \\ 0 & 0 \\ 0 & 0 \end{bmatrix}\,x_1\,x_2 + \begin{bmatrix} 1 & 0 \\ 1 & 1 \\ \hline 3 & 3 \\ 3 & 3 \\ 3 & 3 \end{bmatrix}\,x_2^2 \]
Horizontal concatenation

Consider the following two matrix polynomials:

\[ P_1(\mathbf{x}) = \underbrace{ \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} }_{A_1} + \underbrace{ \begin{bmatrix} 0 & 1 \\ 2 & 3 \end{bmatrix}}_{B_1}\,x_1 + \underbrace{ \begin{bmatrix} 1 & 0 \\ 1 & 1 \end{bmatrix}}_{C_1}\,x_2 % \quad\quad\quad % P_2(\mathbf{x}) = \underbrace{ \begin{bmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \end{bmatrix}}_{A_2}\,x_1 + \underbrace{ \begin{bmatrix} 0 & 1 & 2 \\ 3 & 4 & 5 \end{bmatrix}}_{B_2}\,x_2 \]

Their horizontal concatenation is:

\[ \begin{bmatrix} P_1(\mathbf{x}) & P_2(\mathbf{x}) \end{bmatrix} = \begin{bmatrix} A_1 & \mathbf{0}_{2 \times 3} \end{bmatrix} + \begin{bmatrix} B_1 & A_2 \end{bmatrix}\,x_1 + \begin{bmatrix} C_1 & B_2 \end{bmatrix}\,x_2 \\ \]

which is:

\[ \begin{bmatrix} P_1(\mathbf{x}) & P_2(\mathbf{x}) \end{bmatrix} = \left[\begin{array}{cc|ccc} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \end{array}\right] + \left[\begin{array}{cc|ccc} 0 & 1 & 1 & 1 & 1 \\ 2 & 3 & 1 & 1 & 1 \end{array}\right]\,x_1 + \left[\begin{array}{cc|ccc} 1 & 0 & 0 & 1 & 2 \\ 1 & 1 & 3 & 4 & 5 \end{array}\right]\,x_2 \]

In PolyAny, it's possible to concatenate MatrixPolynomial objects by using the concatenate function from the polyany.functions module.

from polyany.functions import concatenate

Now, you can call concatenate directly.

import polyany.functions as pa

Now, you can call as pa.concatenate.

The concatenate function takes two arguments: the sequence (list or tuple) of matrix polynomials to concatenate, and the axis which informs which type of concatenation (vertical or horizontal) will be performed.

Tip

As a rule of thumb, axis = 0 means a vertical concatenation and axis = 1 a horizontal concatenation.

Warning

To concatenate polynomials, their shapes must be consistent. This means that, in a vertical concatenation (axis = 0), all polynomials must have the same number of columns (dimension 1). When concatenating horizontally (axis = 1), they must have the same number of rows (dimension 0).

Let's reproduce the first example in PolyAny. First, we declare the matrix polynomials:

mpoly1 = MatrixPolynomial(
    [[1, 1], [0, 2]],
    [np.ones((2, 2)), np.tri(2)],
)
mpoly2 = MatrixPolynomial(
    [[0, 0], [0, 2]],
    [np.arange(6).reshape(3, 2), 3 * np.ones((3, 2))],
)

Now, we can concatenate them vertically:

>>> concatenate([mpoly1, mpoly2])  # (1)!
[[0. 0.]    [[1. 1.]            [[1. 0.]
 [0. 0.]     [1. 1.]             [1. 1.]
 [0. 1.]     [0. 0.]             [3. 3.]
 [2. 3.]     [0. 0.]             [3. 3.]
 [4. 5.]] +  [0. 0.]]*x_1*x_2 +  [3. 3.]]*x_2^2
  1. By default, the axis argument is 0.

For the second example:

mpoly1 = MatrixPolynomial(
    [[0, 0], [1, 0], [0, 1]],
    [np.eye(2), np.arange(4).reshape(2, 2), np.tri(2)],
)
mpoly2 = MatrixPolynomial(
    [[1, 0], [0, 1]],
    [np.ones((2, 3)), np.arange(6).reshape(2, 3)],
)

concatenating the polynomials horizontally:

>>> concatenate([mpoly1, mpoly2], axis=1)
[[1. 0. 0. 0. 0.]    [[0. 1. 1. 1. 1.]        [[1. 0. 0. 1. 2.]
 [0. 1. 0. 0. 0.]] +  [2. 3. 1. 1. 1.]]*x_1 +  [1. 1. 3. 4. 5.]]*x_2

Polynomial block

Definition

A polynomial block is formed by concatenating several polynomials both vertically and horizontally into a single polynomial.

Example

Consider the following four polynomials:

\[ P_1(\mathbf{x}) = \underbrace{ \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} }_{A_1}, % \quad\quad\quad % P_2(\mathbf{x}) = \underbrace{ \begin{bmatrix} 0 & 1 & 2 \\ 3 & 4 & 5 \end{bmatrix} }_{A_2} + \underbrace{ \begin{bmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \end{bmatrix} }_{B_2}\,x_1 \]
\[ P_3(\mathbf{x}) = \underbrace{ \begin{bmatrix} 2 & 2 \\ 2 & 2 \\ 2 & 2 \end{bmatrix} }_{A_3}\,x_1, % \quad\quad\quad % P_4(\mathbf{x}) = \underbrace{ \begin{bmatrix} 10 & 11 & 12 \\ 13 & 14 & 15 \\ 16 & 17 & 18 \end{bmatrix} }_{A_4} \]

A possible block of these polynomials is:

\[ \begin{bmatrix} P_1(\mathbf{x}) & P_2(\mathbf{x}) \\ P_3(\mathbf{x}) & P_4(\mathbf{x}) \end{bmatrix} = \begin{bmatrix} A_1 & A_2 \\ \mathbf{0}_{3 \times 2} & A_4 \end{bmatrix} + \begin{bmatrix} \mathbf{0}_{2 \times 2} & B_2 \\ A_3 & \mathbf{0}_{3 \times 3} \end{bmatrix}\,x_1 \]

which is:

\[ \begin{bmatrix} P_1(\mathbf{x}) & P_2(\mathbf{x}) \\ P_3(\mathbf{x}) & P_4(\mathbf{x}) \end{bmatrix} = \left[\begin{array}{cc|ccc} 1 & 0 & 0 & 1 & 2 \\ 0 & 1 & 3 & 4 & 5 \\ \hline 0 & 0 & 10 & 11 & 12 \\ 0 & 0 & 13 & 14 & 15 \\ 0 & 0 & 16 & 17 & 18 \end{array}\right] + \left[\begin{array}{cc|ccc} 0 & 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 & 1 \\ \hline 2 & 2 & 0 & 0 & 0 \\ 2 & 2 & 0 & 0 & 0 \\ 2 & 2 & 0 & 0 & 0 \end{array}\right]\,x_1 \]

To reproduce the example in PolyAny, first we declare the polynomials:

mpoly1 = MatrixPolynomial([[0]], [np.eye(2)])
mpoly2 = MatrixPolynomial([[0], [1]], [np.arange(6).reshape(2, 3), np.ones((2, 3))])
mpoly3 = MatrixPolynomial([[1]], [2 * np.ones((3, 2))])
mpoly4 = MatrixPolynomial([[0]], [np.arange(10, 19).reshape(3, 3)])

Now, we can create the block using the block function from polyany.functions module:

>>> block([[mpoly1, mpoly2], [mpoly3, mpoly4]])
[[ 1.  0.  0.  1.  2.]    [[0. 0. 1. 1. 1.]
 [ 0.  1.  3.  4.  5.]     [0. 0. 1. 1. 1.]
 [ 0.  0. 10. 11. 12.]     [2. 2. 0. 0. 0.]
 [ 0.  0. 13. 14. 15.]     [2. 2. 0. 0. 0.]
 [ 0.  0. 16. 17. 18.]] +  [2. 2. 0. 0. 0.]]*x_1

The block function works by concatenating the inner sequences horizontally and then concatenating the resulting polynomials vertically.

  1. A comparison is made by using np.allclose(), which checks if two arrays are equal within a tolerance. ↩

  2. Python builtins numeric types (int, float) and NumPy scalars. See Scalar ↩↩↩↩↩↩↩↩

  3. Lists, tuples and NumPy 2D-arrays. ↩↩↩