An explicit Wiener-Hopf factorization algorithm for matrix polynomials and its exact realizations within ExactMPF package.

We discuss an explicit algorithm for solving the Wiener-Hopf factorization problem for matrix polynomials. By an exact solution of the problem, we understand the one constructed by a symbolic computation. Since the problem is, generally speaking, unstable, this requirement is crucial to guarantee that the result following from the explicit algorithm is indeed a solution of the original factorization problem. We prove that a matrix polynomial over the field of Gaussian rational numbers admits the exact Wiener-Hopf factorization if and only if its determinant is exactly factorable. Under such a condition, we adapt the explicit algorithm to the exact calculations and develop the ExactMPF package realized within the Maple Software. The package has been extensively tested. Some examples are presented in the paper, while the listing is provided in the electronic supplementary material. If, however, a matrix polynomial does not admit the exact factorization, we clarify a notion of the numerical (or approximate) factorization that can be constructed by following the explicit factorization algorithm. We highlight possible obstacles on the way and discuss a level of confidence in the final result in the case of an unstable set of partial indices. The full listing of the package ExactMPF is given in the electronic supplementary material.

Introduction and motivation

The factorization problem has proven to be an extremely efficient tool in mathematics (for example, in the theory of linear operators [1,-3]) and instrumental in various branches of physics, mechanics and engineering (see, for example, [4,-8]).

Unfortunately, in contrast with the scalar problems, factorization of arbitrary matrix functions remains an unresolved problem [9,10] (especially, when one speaks of its effective applications). We start with definitions, preliminary results and than formulate the problem considered in this paper.

Let be a simple piecewise smooth closed contour in the complex plane bounding the domain . We can assume that . The complement of in the closed complex plane, , is denoted by . Let be matrix function on . In what follows we assume that is continuous and invertible on .

The representation of a matrix function is called the right Wiener–Hopf factorization of the matrix function [1]. Here, are continuous and invertible matrix functions on that admit analytic continuation into and their continuations are invertible into the respective domains; matrix function , where integers are called the right partial indices of . It is always possible to assume that . Moreover, , where is the winding number (or the Cauchy index) of the determinant of the matrix .

Similarly, the left Wiener–Hopf factorization of is the representation [1] where have the same properties as , , and are called the left partial indices of the matrix function (). Note that, if the factorizations exist, their left and right indices are usually different sets of integers and constructions of the right and left factorizations are usually considered as two separate problems. On the contrary, a method proposed in this work requires simultaneous considerations of both factorizations.

In general, a matrix function with continuous entries does not admit the Wiener–Hopf factorizations with factors in the same class. It is known that if the elements of satisfy the Hölder condition on or belong to the Wiener algebra , where is the unit circle in the complex plane, then the Wiener–Hopf factorizations exist (see [1]).

The scalar factorization problem for a non-singular function can always be solved explicitly, see, for example, [11] and [12]. Unfortunately, explicit formulae for the factors , and the partial indices of an arbitrary matrix function do not exist at present.

Fortunately, there are classes of matrix functions for which an explicit (or constructive) solution to factorize the problem have been found. By this, we understand that there are explicit formulae allowing us to formally construct the partial indices and the factors of the matrix function . For example, the problem for matrix polynomials has been solved explicitly [-1,13,15]. For meromorphic matrix functions, there also exits an explicit solution [13]. A detailed review of constructive methods for the factorization problem is presented in the work [10] (see, also [16]).

Unfortunately, the factorization problem for matrix functions is generally speaking unstable. More precisely, if the difference between the larger and the smaller partial indices is greater than one, then the set of the partial indices can change with an arbitrary small perturbation of the initial matrix [12,17]. Moreover, if left (right) factorization appears to be stable, this does not mean that the right (left) preserves that property [1]. There is also a technical difficulty embedded into the problem as the factors , are not found uniquely, while how to properly normalize the algorithm to guarantee the uniqueness in the chosen class is still open. On the other hand (see, for example, [18], theorem 6.15), if two matrix functions and are close enough, and their partial indices coincide, then it is always possible to find a factorization of for which factors are close enough to the factors of . Recently, a few constructive approximate techniques for factorization of specific matrix functions have been proposed [19,-22], under the condition that their sets of partial indices are known in advance. Accompanied stability analysis has also been performed in the case of unstable partial indices.

Therefore, the key step in construction of any (exact or approximate) factorization of a given matrix function (if the factorization exists) is to provide a procedure/formula allowing for calculation of the partial indices. However, even having an explicit route to achieve the goal, it often requires a numerical realization. And here, one may meet with a hidden unexpected obstacle: the procedure can again be unstable. Below we discuss two simple examples highlighting the issue.

Consider two matrix polynomials (the only elements are slightly different there): To factorize explicitly the polynomials by the method developed in [14], we need first to factorize their determinants. For the first matrix function, we have and the factorization of this polynomial can only be done approximately. Hence, all subsequent steps for the explicit factorization may be unstable. To apply the effective criterion of stability for partial indices [23], we have to verify that a certain matrix with float point entries, built using the matrix , has the full rank. In turn, this is an ill-condition problem itself if the condition number of the matrix is large enough ([24], ch. 2). Thus, within the framework of this particular and well-developed explicit method, we cannot always even determine for sure whether the partial indices are stable or not. All other methods for matrix polynomials known in the literature [15] have the same obstacle when used to factorize the matrix .

A completely different situation arises with the polynomial . Now the determinant also has its roots that can be found only numerically (and thus, approximately). However, it can be presented as a product of irreducible polynomials and all roots of lie in an interior of the unit circle, while all roots of the polynomial lie in the domain exterior. We can verify these statements without calculating the roots themselves by using Schur’s test (see §3). Hence, is the Wiener–Hopf factorization of the determinant, , over the field of the Gaussian rational numbers , since the coefficients of the matrix polynomial belong to this field. In what follows, we show that all steps of the explicit factorization algorithm can be performed in the exact arithmetic and, as a result, the instability issue does not arise in this case.

These examples demonstrate that By the explicit (or constructive) solution of the factorization problem, we understand a clearly defined algorithmic procedure that should definitely terminate after a finite number of steps. However, when one implements a specific explicit factorization algorithm, each consequent step can be executed exactly or approximately (numerically).

existence of an explicit procedure for a particular matrix factorization does not guarantee that it can be successfully used in practice;

knowledge of exact values of zeros of the polynomial is not an essential feature (useful though) for execution of the exact factorization;

the notions of explicit, exact and numerical factorizations should be clearly distinguished.

To avoid the difficulties associated with possible instability, one can chose the input data belonging to the field and perform all steps of the explicit algorithm in the exact arithmetic. In such a case, we find that there exists an exact solution of the factorization problem. Unfortunately, as has been discussed above, this resolution is not always available even when in possession of a few different explicit algorithms.

If a factorization problem does not admit the exact solution, the only choice that remains is to construct a numerical (approximate) factorization. The latter, in view of possible instability, may represent a real challenge for execution and all computations should be done with special care and supplemented by accurate analysis of the partial indices preservation and a numerical stability.

Thus, the key question arising during the factorization of a matrix function is: how to compute the partial indices with a full confidence?

There are classes of matrix functions for which the partial indices are known in advance. For example, positive definite matrix functions [12] have all their partial indices equal to zero, while, for matrix functions not containing zero in their numerical ranges [25], they are equal (but not necessarily zeros). Such sets of the partial indices are stable with respect to small perturbation of the matrix functions [12], and, as a result, the remaining task in the factorization process (computation of the factors) is less problematic and any appropriately chosen numerical algorithm can be performed and justified by standard techniques.

In this study, we discuss factorization of matrix polynomials on the unit circle with particular attention to all these details. The structure of the paper is as follows. In §1, we describe a well-known explicit method for constructing a factorization of a matrix polynomial [14]. We formulated this method algorithmically in order to indicate explicitly all the steps that must be taken to construct the factorization. This will be convenient for further attempts to implement our method numerically. However, in this section, we do not consider computational problems related to our method. Therefore, no analysis of the accuracy (forward/backward error) of the explicit algorithm is carried out in this section. The basic tool for this factorization is a method of essential polynomials [26].

In §3, we show that all steps of the explicit algorithm for factorization of a matrix polynomial, , over the Gaussian rational field is implemented in the exact arithmetic if and only if its determinant, , allows exact factorization in this field. The solution of the latter is presented in theorem 3.1. We used representation a scalar polynomial as the product of its irreducible factors and the Schur algorithm to test whether each factor has its roots siting inside (or outside) the unit circle. Further, all steps of the explicit matrix factorization can be performed in the exact arithmetic. The algorithm for exact evaluation of the essential polynomials is a foundation stone of theorem 3.2 which is the main result of this section. The formal description of this algorithm is rather complicated and will be presented in the separate work.

In §4, we implemented the exact factorization algorithm as the package ExactMPF in the computer algebra system Maple.[1] This package consists of the following basic procedures: SchurTest, SolverExactPF, ExactMLC, seqTk, NumRank, ExactFEP. The main procedure is SolverExactMP. The purposes of these procedures are shortly described in the same section, while the full listing of the package ExactMPF is given in the electronic supplementary material. Several examples of the exact factorizations performed with use of the package ExactMPF are presented in §5 with the aim to demonstrate the package abilities and to study peculiarities of the factorization for special classes of matrix functions. Other tests are used to analyse stability of the computations (see examples 5.1 and 5.2). All those and other examples are included in the electronic supplementary material with more details.

In the general case, the explicit algorithm can only be realized approximately. This poses a real challenge for users due to the instability of the problem. Our numerical experiments with ExactMPF have demonstrated that even in the case of unstable sets of partial indices there is a good chance to obtain an approximate factorization. Namely, to accelerate the work of ExactMPF, we provide an option ‘symnum’ into the SolverExactMP. If this option is turn on, the package calculates the ranks of matrices by the singular value decomposition (SVD) method [24,27]. These numerical ranks are used to find partial indices of the initial matrix polynomial. Further calculations are performed in the exact arithmetic. This allows for a verification of the partial indices finding by SVD versus the exact ones: (are the sets of the partial indices the same or not?). It turned out that for many examples, even with unstable sets of the partial indices, such a numerical approach leads to correct results. The partial indices found with the SVD method are called partial -indices of the polynomial. However, the question whether or when they do coincide remains open. The opportunity to use SVD for numerical calculations of partial indices in a more general context is discussed in §6. Finally, in the Conclusion, we underline the advances made in this paper and discuss further possible developments and applications.

An explicit Wiener–Hopf factorization for matrix polynomials

In this section, following the results of the work [14], we present an explicit algorithm for the factorization of matrix polynomials. In accordance with our understanding of the notion of explicit solution (see the Introduction), we describe below those finite numbers of steps that, after their execution, allows one to construct a factorization. We focus on the algorithmic structure of this explicit procedure that should be implemented within a respective numerical code. We would like to emphasize here that this section is not devoted to any study of possible computational problems that may appear during the algorithm implementation into a specific environment.

We start with Input data, that is, a matrix polynomial of a degree , , represented by the set of its consecutive coefficients, , , and defines the matrix size. We assume that the matrix polynomial, , is invertible on the unit circle .

Construction of the Wiener–Hopf factorization of the scalar polynomial .

We denote . Since , the determinant admits the Wiener–Hopf factorization relative to , that is unique with the additional condition at infinity for the polynomial . In the sequel, we use, in fact, only one of the factors, namely,

It should be noted that any explicit factorization algorithm will use the calculation of the determinant of the factorizable matrix function. The authors realize that the complexity for computing the determinant of a matrix of order is . Fortunately, in applications this happens for matrix function of a small order. Moreover, in §3, we use the matrix determinant in conditions of the exact calculations.

Here, is the Cauchy index of . There is Gakhov’s formula for this, which is convenient for numerical calculations. However, in our case, we can find this index by the exact calculations (see theorem 3.1).

Evaluation of a finite set of Laurent coefficients of the matrix function

Laurent coefficients of the auxiliary rational matrix function are computed in terms of matrix coefficients of the original matrix polynomial, , and the coefficients of the scalar polynomial, , defined on the previous step. The rational function is analytic in a domain for some . Let us expand and in the Laurent series in this domain. Coefficients and , in the series above, can be found from the following relations: and where the factor has been computed on the previous step in (2.2). In what follows, we need only a finite number of the coefficients, , for .

Let be the finite sequence of complex matrices which are the Laurent coefficients of . We form a finite family of the block Toeplitz matrices of finite sizes The next step is a corner stone in implementation of the explicit algorithm proposed in [14].

Finding the ranks and the null spaces for of block Toeplitz matrices .

Let () be the right (left) kernels of the matrices . Recall that the right and left kernels of a complex block matrix consisting of blocks are defined by the relations .

Let us represent an element in a block form , and define a column-valued generating polynomial in related to as

By , , we denote the space of generating vector polynomials for vectors in . Put and let be –dimensional space of all column-valued polynomials whose degrees are not greater than .

Repeating the same line of reasoning, we define the space , of the row-valued generating polynomials in for the rows from .

As a result of this step, the set of ranks of the matrices defined in (2.6), and bases of the spaces , , , has been computed. All those elements, values and sets are of finite dimensions.

In general, computing rank and calculation of basis elements of the null-spaces are unstable processes. In our work, we will find the rank and null-space of a matrix with the help of calculation in the exact arithmetic.

Computation of indices and factorization essential polynomials.

The main tools for computations of the partial indices and factors in the factorization of matrix functions are the so-called indices and essential polynomials of the sequence (see [14,26]). Below we give necessary definitions and discuss how to compute them.

By , we denote a dimension of the right kernel and introduce the following integers: for . A sequence is called regular if and . For a regular sequence, we have (see [14,26]) Since a monotone integer sequence is piecewise constant, then there are integers such that The absence of the th row here means that .

Definition 2.1.

The integers defined by the relations (2.8) are the indices of the sequence .

It turns out that (see [14,26]) Furthermore, we define the right essential polynomials of the sequence . Note that and are subspaces of as it follows from the definition of the spaces . The dimension, , of the complement of in is equal to .

Then, equations (2.8) imply that if and only if , . Moreover, in this case, is equal to the multiplicity, , of the index . Therefore, for , we have and for

Definition 2.2.

Any polynomials forming a basis for a complement are called right essential polynomials of the sequence corresponding to the index .

As a result, we have defined indices and right essential polynomials for any regular sequence . Similarly, we can define the left essential polynomials of the sequence .

By theorem 3.1 from [14], the sequence defined in equation (2.5) is regular and there exist respective essential polynomials ; such that

(i) the constant terms of the polynomials are equal to zero,

(ii) the leading terms of the polynomials are equal to zero.

Definition 2.3.

The essential polynomials satisfying the conditions (i)–(ii) are called the factorization essential polynomials of the sequence.

Step 5. Construction of the factorization of the initial matrix polynomial,

Now we are in a position to deliver a final result on the explicit Wiener–Hopf factorization of a matrix polynomial (see theorem 3.2 in [14]). It is crucial to note here that, in contrast to any other methods, if it exists for a particular matrix function, our approach requires simultaneous consideration of both the right and the left factorizations. In the statement of this theorem, we have corrected the misprints appearing in the formulae for the factors , in theorem 3.2 of [14].

Theorem 2.4.

Let be the indices and are the right (left) factorization essential polynomials of the sequence . Let us introduce the matrix functions and Then the left () and right () partial indices and the factors (, ) of the respective factorizations of the matrix polynomial are defined by the formulae

It follows from this theorem that the indices of the sequence , in addition to (2.9), satisfy conditions

Those five steps above represent an explicit algorithm, since their implementation requires only using the linear algebra tools for objects of finite dimensions.

Remark 2.5.

The main difficulty in the practical implementation of the aforementioned formulae is that the sequence can be found, in general cases, only numerically (and thus, the result is computed with some accuracy). Therefore, in such cases, the algorithm would produce the factorization of a matrix polynomial close enough to the original matrix polynomial with all consequences of such replacement (possible instabilities as discussed above).

Remark 2.6.

If, by any means, one believes that the partial indices of a matrix polynomial in question are stable, some computations from Step 4 can be avoided (see example in §6).

An exact realization of the explicit algorithm

The analysis of the algorithm, carried out in the previous section, shows the reason for the instability of the factorization problem for matrix polynomials. The explanation is in the instability of the problems of finding the rank and null-space of a matrix. Since finding the rank is the key problem for our algorithm (and also for the Gohberg–Lerer–Rodman algorithm [15]), in the general case, an explicit algorithm cannot be implemented numerically.

However, if Steps (1)–(5) can be carried out in exact arithmetic, i.e. the error-free calculations are used, problems with the factorization instability, with the accuracy of calculation and rounding errors do not arise. It should be noted that symbolic calculations have the significant drawback that they often require significant computational resources. In our case, an increase in the size of the matrix function and the degree of the matrix polynomial can lead to a significant slowdown in the algorithm. Fortunately, in most of the applications, matrix functions of the second or third order mainly appear, and the only reason for a possible slowdown is a large degree of the matrix polynomial. In table 1, we illustrate this feature in detail.

Table 1

Execution time, , (in seconds) for factorization of the matrix polynomial , defined in (5.4). The letter ‘F’ indicates a failure occurred during the execution, while the symbol ‘–’ highlights that the computations, performed in the symbolic mode, were unable to complete the task within 1 h.

						t_k (symnum)
k	n_k	p_k	ϰ_k	ϰ_2 (T_0)	t_k (sym)	digits=5	digits=10	digits=15
1	4	6	3	1.10×10	0.28	0.4	0.31	3.44
2	8	12	6	7.95×10	0.81	1.0	0.61	0.62
3	12	18	9	6.04×10^2	2.62	1.4	1.77	1.80
4	16	24	12	5.78×10^3	7.19	F	2.94	2.95
5	20	30	15	5.90×10^4	23.78	F	7.95	8.83
6	24	36	18	5.95×10^5	55.11	F	10.66	13.26
7	28	42	21	5.66×10^6	147.23	F	35.80	38.11
8	32	48	24	5.13×10^7	288.67	F	44.89	42.19
9	36	54	27	4.65×10^8	644.03	F	F	135.55
10	40	60	30	4.44×10^9	1103.06	F	F	150.59
11	44	66	33	4.53×10^10	2204.62	F	F	453.78
12	48	72	36	4.66×10^11	—	F	F	444.05
13	52	78	39	4.57×10^12	—	F	F	F

The main result is formulated in theorem 3.2. It gives a necessary and sufficient condition for the exact factorization of a matrix polynomial , relying on the same property of its determinant, . Below, we discuss a practical implementation of those five steps of the explicit algorithm in this particular arithmetic.

Input data. A matrix polynomial , , of degree that is invertible on the unit circle . The entries of the matrix coefficients belong to the field .

Step 1. In order to construct the Wiener–Hopf factorization of the scalar polynomial in exact arithmetic, it is sufficient to require that all roots of belong to the field . However, this condition is not necessary. There exists a weaker necessary and sufficient condition allowing exact Wiener–Hopf factorization of over the field . The idea here is to factorize the polynomial as a product of irreducible factors, testing (for example, with help of the Schur algorithm) whether the roots of each factor lie in the interior or exterior of the unit disk.

First, we give the criterion in terms of coefficients of the irreducible factors of the polynomial. As has been mentioned above, we assume that no zero lies on the unit circle. Note that, in the opposite case, the Schur test fails during the execution stage. This fact explains why we have not concentrated on how to verify the condition on the first step of the algorithm.

Let be a polynomial with complex coefficients, . We denote , where is the complex conjugate of .

For the given , we define the triangular Toeplitz matrices and form the Hermitian matrix . Here, () is Hermitian conjugate of ().

Let be a polynomial over the field . It can be represented as a product of irreducible monic polynomials [28]:

Theorem 3.1.

The polynomial admits the exact Wiener–Hopf factorization if and only if, for every polynomial , , one of the matrices or is positive definite.

Suppose these conditions are fulfilled. Denote by the product of those irreducible factors in decomposition (3.1) for which the matrices are positive definite. Then, the polynomial is the product of the remaining irreducible factors and . Let us define: , and . Then the exact Wiener–Hopf factorization of the polynomial is

Proof.

It is clear that the polynomial admits the exact Wiener–Hopf factorization if any of the irreducible factors has all its roots lying inside or outside the unit circle. By the criterion Schur–Cohn [29,30], the polynomial has all its roots lying inside the unit circle if and only if the matrix is positive definite. Similarly, has all roots lying outside the unit circle if has all its roots lying inside the unit circle, that is is positive definite. The representation (3.2) is the Wiener–Hopf factorization that is normalized by the condition .

For practical computation, it is more convenient to use the Schur test [31] instead of checking the positive definiteness of matrices. This test states (see, for example, [29]) that

a polynomial

, and

the polynomial determined from the equality possesses the same property.

Since , then in a finite number of steps, we arrive at a polynomial of degree 0 or find that the polynomial does not have the required property.

In order to verify that all roots of lie outside the unit circle, it is necessary to apply the Schur test to the polynomial . Note that the above algorithm automatically verifies the invertibility of the matrix polynomial on the unit circle .

Step 2. Since the polynomial admits the exact factorization the coefficients of the polynomial in belong to the field . The function is analytic in the domain , for some . Its Laurent coefficients are obtained by equation (2.4) in exact computations. Since only a finite number () of the coefficients are needed and all operations are arithmetic, Step 2 is therefore performed in the exact manner.

Step 3. The matrices defined in equation (2.5) have the entries belonging to , i.e. the sequence is exactly computed. In Step 3 of the explicit algorithm, we find ranks and null spaces of block Toeplitz matrices for . This step can be realized by Gaussian elimination. Since entries of the matrices belong to , the Gaussian elimination can be implemented in exact arithmetic [28].

Step 4. Obtaining the indices and the factorization essential polynomials of the sequence , defined in (2.8) and (2.11), is a foundation stone of the method. Since the ranks of the matrices have been calculated exactly, the differences, , and hence the indices, , are also found exactly.

To determine the essential polynomials, we need bases of the complements , , (see equation (2.10) and definition 2.2) and the factorization essential polynomials (see definition 2.3). It is intuitively clear that these polynomials can be found exactly using the methods of linear algebra. However, a rigorous proof of this fact is rather complicated and requires a careful study of all steps of the algorithm for determination of the factorization essential polynomials. It will be presented in the separate work.

Step 5. The main result of this section is the following.

Theorem 3.2.

Let be a matrix polynomial over the field of Gaussian rationales . Then admits exact Wiener–Hopf factorization if the polynomial is exactly factorable over .

If admits the exact left Wiener–Hopf factorization , then is the exact factorization of .

Conversely, if the polynomial admits the exact factorization then all preliminary Steps 1–4 can be carried out exactly. It remains to use the formulae (2.12)–(2.14).

Implementation of the exact algorithm in Maple’s code: ExactMPF package

The exact factorization algorithm has been realized in the computer algebra system Maple using its linear algebra utilities as the package ExactMPF. The package can work in the version Maple 11 or higher. Clearly, one can do this in any other suitable software (for example, Mathematica[2] ). The main procedure of the package is SolverExactMPF. During the execution, it calls a number of subroutines. We list the basic ExactMPF’s procedures:

SchurTest verifies whether all polynomial roots lie in the unit circle (see §3, Step 1).

SolverExactPF verifies whether a scalar polynomial admits the exact Wiener–Hopf factorization and builds it (see §3, Step 1).

ExactMLC computes the matrix Laurent coefficients for the rational matrix function (see equation (2.5) and §3, Step 2).

seqTk builds the sequence of the matrices , (§3, Step 3).

ExactFEP finds indices and factorization essential polynomials of the sequence .

The main procedure SolverExactMPF finds the exact factorizations of a Laurent matrix polynomial . It also checks and validates the input entered into the procedure. The following validations are performed: If the validation fails, the program is interrupted.

are elements of the matrix polynomials in , ?

are these polynomials over the field ?

does admit the exact Wiener–Hopf factorization?

To access ExactMPF use the commands.

> read("ExactMPF.txt");

> with(ExactMPF);

> with(LinearAlgebra);

To obtain the factorizations of we run the procedure SolverExactMPF with argument .

lplus, dl, lminus, rminus, dr, rplus := SolverExactMPF(a):

The procedure SolverExactMPF returns the factors lplus, dl, lminus of the left factorization and the factors rminus, dr, rplus of the right factorization.

This procedure also automatically performs a few tests for the constructed factorization

it finds the sum of the partial indices,

checks whether the factors () are matrix polynomials in (in ) and

whether , are constants.

A very peculiar option symnum is also built into the procedure SolverExactMPF. It can be evaluated by the command:

lplus, dl, lminus, rminus, dr, rplus := SolverExactMPF(a,"symnum"):

In this mode, SolverExactMPF converts the matrices to matrices with floating point entries to the precision given by the global variable digits. Then, SolverExactMPF uses the procedure NumRank included in the package ExactMPF.

NumRank finds a numerical rank of a matrix with float point entries. A good practice is to use a numerical rank, which is defined with the help of the SVD method [24,27]. Recall some necessary definitions.

Let be a matrix and , , its singular numbers. In the case of exact computation, if , then and . The number is called the spectral condition number of the rectangular matrix [27,32].

It is known that perturbations of the singular values are of the same order of magnitude as the perturbations of the matrix entries. Thus, computation of the singular values is a stable procedure and this property has led to the definition of a numerical -rank (see [24,27]).

A numerical -rank of the matrix function is the number of singular values greater than some tolerance , i.e. the singular values that are less than or equal to considered as perturbations of the zero singular values.

Thus, the tolerance specifies a way of separating the non-zero singular numbers from those assumed to be zero before the perturbation. Hence, the notion of -rank is useful when there exists a well-defined gap between non-zero and zero singular values of the matrix. Even for a matrix of the full rank, this requirement may not be fulfilled (if the matrix is ill-conditioned related to the spectral condition number). Choice of the parameter may be very sensitive and thus it is the most delicate and difficult task to decide upon in calculation of the numerical -rank. In the procedure NumRank, a default value of the tolerance is .

Further, -ranks of the block Toeplitz matrices are used to define -indices of the sequence instead of the exact indices of this sequence. All remaining calculations are performed in the exact arithmetic. If the output is validated successfully, the replacement of the indices by -indices may be considered as justified. Using this mode, one can significantly speed up the calculations. Moreover, this mode allows us to perform numerical experiments to test a coincidence of the partial indices and the partial -indices that were found by -indices of the sequence . In the next section, we perform such experiments.

After this short description of the ExactMPF package, we come back to matrix polynomials and discussed in the Introduction and try to factorize them exactly using the package.

Example 4.1.

Consider first the matrix function defined in (1.2) and execute

> SolverExactMPF(a):

‘Failure: The exact polynomial factorization doesn’t exist’

‘The factorization process is interrupted’

Example 4.2.

Consider the second matrix polynomial from (1.2) that allows for the exact factorization of its determinant over the field . The procedure SolverExactMPF gives in this case the following expressions for the factors of :

> lplus; dl; lminus;

> rminus; dr; rplus;

The execution time is 0.265 s.

Note that the left and right partial indices do not coincide. Moreover, the set of partial indices in the left factorization is stable, while that for the right factorization is unstable. In the next section, we carry out more extensive numerical experiments with ExactMPF.

Examples of matrix factorization using the ExactMPF package

All computations have been performed on a home desktop computer HP with Intel(R) Core(TM)i3-415T CPU, 3.00 GHz, 4G RAM, operating the system Windows 10.

The package ExactMPF can be used for study of the explicit factorization of matrix polynomials. In the electronic supplementary material, we include numerous tests performed with use of the ExactMPF package, containing solutions of all examples from the main body of the paper.

Factorization of triangular matrices

First, we consider factorization of triangular matrix functions. In [33], the right factorization of triangular matrix functions is studied in detail and the stability analysis of this problem is carried out. It has been proven that, for a non-trivial case, the factorization of is reduced to the factorization of the following Laurent matrix polynomial where . If , factorization of such a matrix can be naturally constructed with the help of our package.

Example 5.1.

Consider Then, after execution of SolverExactMPF, we get

> rminus, dr, rplus;

> lplus, dl, lminus;

The execution time is 0.547 s. We note that in this case, the left factorization can be easily constructed without the package, but with another order for the left partial indices

As is clear from the presented example, the indices of the diagonal elements in the triangular matrix do not always coincide with its partial indices [33].

The factorization of triangular matrix functions of an arbitrary order requires considerations of a number of special cases, depending on the relation between the indices of the diagonal elements [34]. The procedure SolverExactMPF can be useful for experimenting with such matrices.

Example 5.2.

Let

Here, the indices of diagonal elements are , , . SolverExactMPF returns the following data:

> lplus, dl, lminus;

Since the expression for the right factors are cumbersome (see the electronic supplementary material for full details), we present only the middle factor

> dr

The executing time is 0.563 s.

Again, both the left and right partial indices of this matrix function do not coincide with the indices of the diagonal elements.

Factorization of a matrix polynomial with complex coefficients

In the following example, we construct the exact factorization of a matrix polynomial with complex coefficients.

Example 5.3.

Then, we get

> lplus, dl, lminus, dr;

The execution time is 0.438 s.

For the same reason as above, we present the factors only for the left factorization.

Experiments with the -factorization of matrix polynomials

As have been already underlined, the main unsolved problem in the approximate factorization of matrix polynomials is a verification whether its partial -indices coincide with the partial indices of the original matrix (both left and right ones). Since it has been impossible so far to prove any definite result theoretically, we carry out here numerical experiments using our designed solver SolverExactMPF in its symbolic-numeric mode symnum for the cases when exact factorization can be performed.

In the next example, we introduce a sequence of exact factorizable matrix polynomials with increasing degrees. We use the full capability of SolverExactMPF in both its modes (exact and symbolic-numeric) to calculate their partial indices and the partial -indices (see §4). Then we carry out all subsequent calculations in the exact arithmetic. If the obtained factors , and the partial indices pass all the necessary verifications, then the factorizations are constructed correctly and the replacements of the ranks by their numerical -ranks is justified. Moreover, using option symnum mode gives an additional bonus: it significantly speeds up the computations.

The purpose of our experiments is to test applicability of the notion of the -rank for the factorization of matrix polynomials. A priori, it is clear that the method is ineffective if the task of computation of the -ranks for the sequence of block Toeplitz matrices is extremely ill-conditioned. Let us define an analogue of the condition number for this task.

As usual, if is a rectangular matrix, then , , is the condition number of with respect to the spectral matrix norm. We say that is the condition number of the matrix sequence .

Since the -rank can be found only numerically, we are able only to estimate the value of . It is easy to verify that . Hence, we can explicitly calculate and can find a lower bound for . As has been mentioned above, a choice of parameter , to be used in -factorization, is not trivial. We use in our analysis the value as it is the same as the default tolerance in Matlab’s rank command. We can prove that where is the Frobenius norm of a matrix.

Taking into account (5.2), we see that chosen as the tolerance is controlled by the value of parameter digits.

In the next example, we will see how the correctness of such an approach depends on the choice of the parameter digits.

Example 5.4.

The matrix polynomial admits the exact factorization. Let us consider matrix polynomials and denote , , . We also record the execution time, , required by the package to compute factorization of the matrix polynomial . Three different machine precisions: are used for the testing.

The results of the numerical experiment are shown in table 1. Expressions for the factors , are very cumbersome and are omitted.

The test has demonstrated that the symbolic part of the ExactMPF package is able to deliver the results in reasonable time for and the matrix polynomials have the stable sets of their partial left and right indices. Those are From table 1, we can also see that the executions in the sym-mode, when the ranks are calculated in exact arithmetic, become very slow for large values of and the use of symbolic-numeric calculation (symnum-mode) can greatly accelerate the computations.

In the cases indicated by the letter F in table 1, the procedure NumRank has found the ranks of incorrectly (the validation of the sum of the indices for the sequence fails), and the ExactMPF package interrupts computations because of the failure. The reason for this is that the condition number was too large. In this case, increasing parameter digits (that decreases the threshold for numerical tolerance, ) makes it possible to continue correct calculations. Note that the partial indices computed for with use of the -factorization follow the same rule as in (5.5). Unfortunately, this fact does not formally allow us to conclude (with 100% confidence) that the exact partial indices would be the same.

The next example shows that symnum mode is able to correctly calculate the partial indices even in unstable cases.

Example 5.5.

Let us consider a matrix polynomial The exact computations using SolverExactMPF give and , . Moreover, in this case, .

Next, let us construct another matrix polynomial . For this one, SolverExactMPF returns the following values of the partial indices: , and , . In the symnum mode with the parameter digits=10, the package was unable to find the correct values since the tests on the sums of indicies are not successful and the factorization process is interrupted. In this case, and, as a result, the condition number is too large. If, however, we set digits = 15, the package returns, in the symnum mode, the correct values of the partial indices.

Thus, our numerical experiments with several matrix polynomials (see electronic supplementary material) have shown that, for the case of not too large condition numbers, calculations with use of both the symbolic and numeric options of the ExactMPF package give identical partial indices and partial -indices, even in the unstable cases. This fact gives hope for a stable construction of an approximate factorization but requires a solid theoretical proof.

On numerical realization of the explicit algorithm

In previous sections, we considered the problem of the Wiener–Hopf factorization for matrix polynomial over the field . In this case, As a result, we have constructed the ExactMPF package for the exact solution of the factorization problem, that returns the solution under condition that the determinate of the matrix polynomial is exactly factorizable or, in the opposite case, interrupts the computations. In this case, the stability and accuracy problems do not arise. However, a natural question arises: can we deliver an approximate (numerical) factorization of such a matrix? If the conditions for the existence of the exact factorization are not fulfilled, then it is not yet possible to build any algorithm for the factorization problem that satisfies the entire canon of computational mathematics (with carrying stability analysis, obtaining forward/backward errors and accuracy for approximate factorization).

the input data are given with infinite precision, and

there is a possibility to use the exact (error-free) calculations.

Here, we highlight possible obstacles when implementing the explicit algorithm (Steps 1–5) numerically without relating the steps to the ExactMPF package, in the case, when, at least, one of the conditions (i), (ii) is not fulfilled.

Input data. We consider a matrix polynomial of the degree and size that is invertible on the unit circle . The entries of the matrix coefficients, , belong to the field . They can be represented by floating point numbers with the accuracy (finite precision) or the entries belong to the field (infinite precision).

Step 1. Since one of the conditions (i), (ii) is not fulfilled, we can solve the factorization problem for the function only approximately. If the input data are given with infinite precision we convert them in the floating point format with a desired accuracy .

A standard way for the factorization of a scalar polynomial, recommended by all explicit methods, suggests using approximate values of the polynomial roots. However, those roots are, in general, not well-conditioned functions of polynomial coefficients (see, e.g. [35]), and vice versa, the coefficients of a polynomial are also not well-conditioned functions of their roots [36]. The latter means that searching numerically for zeros of the polynomial to construct its factorization (2.1) may lead to the algorithm instability.

In order to solve the factorization problem for scalar polynomials over the field , a dedicated package PolynomialFactorization (in the Maple Software) has been developed by Adukov [37]. It factorizes an input polynomial based on the indices and essential polynomials technique. Input data for its main procedure SolverPF is a scalar polynomial and the accuracy of the polynomial coefficients. SolverPF returns factors , of the Wiener–Hopf factorization of and their guaranteed accuracy . The accuracy, , calculated by the procedure SolverPF depends on a type of polynomial and the accuracy of the input data.

Note, that the PolynomialFactorization package is not included in the ExactMPF package, and is different from the standard ‘Factoring a Polynomial’ Maple command Factor.[3] Its main feature is: it does not rely on computation of the polynomial’s roots.

Step 2. Now we can calculate the Laurent coefficients of using given matrix coefficients of and taking into account their accuracy.

Hence, we can find the sequence with some accuracy as described in §2. This step does not contain any technical difficulties. Note that, if the input data are given with infinite precision, i.e. belong to the field , then we can calculate with any reasonable accuracy. Otherwise, the accuracy is restricted by finite accuracy of input data.

Step 3. This is a crucial point in the factorization problem. In the exact arithmetic, the respective calculations do not meet any problem (except of computational time). However, if the entries of the matrices are known only approximately, evaluations of the ranks and null spaces of might be unstable. This is a reason for the instability of the factorization problem itself.

In §4, we have discussed the procedure NumRank from ExactMPF in detail. It uses the SVD method for finding a numerical rank of matrices with floating point entries. The procedure can also be used in this case.

For definition of the essential polynomials, we need a method for evaluation of a basis of the null space for matrices with floating point entries. This can also be done with the SVD method [24,27]. Computation practice shows that this method is the most reliable way to numerically find the rank and the basis for the null space of a matrix. Here, however, the main problem is how to select the computational tolerance . Unfortunately, in general, there is no rigorous way to do this. As an initial choice, we can recommend using , if one works in the Maple environment.

Step 4. Since entries of the matrices () are, in general, floating point numbers, we need to compute the numerical -rank of the sequence instead of exact. The indices of the sequence found by -ranks are respectively called -indices.

Next, to construct an approximate Wiener–Hopf factorization of the matrix polynomial, we need to find the factorization essential polynomial of the sequence . In §3, we describe an algorithm for construction of these polynomials and prove that, for input data with infinite precision, the algorithm can be realized in the exact arithmetic. An open problem, however, is how to adapt the algorithm for input data with a finite precision. In such cases, we may use a basis of null space for matrix () via the SVD method with the tolerance (-bases). The factorization essential polynomials, constructed with the help of -bases, are called factorization essential -polynomials.

Step 5. For the fixed tolerance , the integers , , computed by formulae (2.12) with the help of -indices of the sequence are called partial -indices of the matrix polynomial . The factorization of constructed by formulae (2.12)–(2.14) with the help the partial -indices and the factorization essential -polynomials will be called -factorization of the matrix polynomial .

The key question here remains open:

Does the

To solve this problem, we must find answers to the following basic questions:

does a tolerance exist such that the partial -indices of coincide with the matrix partial indices? How can this tolerance be found?

can the algorithm for constructing factorization essential polynomials from §3 be adapted to find the factorization essential -polynomials?

In order to carry out numerical experiments, we should develop a new appropriate software. At the present time, we have in our hands the package PolynomialFactorization from [37] and the procedure NumRank from this ExactMPF package. Since the procedure for finding the factorization essential -polynomials is still unavailable, in the following example, we restrict ourselves to the evaluation of partial -indices only.

We consider the matrix polynomial from the Introduction in (1.2). Note that the input data have infinite precision, while the determinant can be factorized only approximately. Let us first construct a factorization of with the help of the package PolynomialFactorization to obtain and with a precision . This precision was returned by the package PolynomialFactorization proceeded under the Maple value of and accuracy .

Thus, an approximate factorization (2.1) has been constructed with the index and the factors and

The Laurent coefficients of are found directly from the recurrence formulae (2.4) Calculation of the matrix Laurent coefficients by formulae (2.5) gives the following result: and

In order to verify the stability criterion for matrix polynomials (theorem 3.1 in [23]), we need to compute the -ranks of matrices

This can be done with use of the function SingularValues(A) of the MapleSoftware that returns the singular values of a matrix . For the matrix , it gives the following singular values: and, therefore, for .

Since the gap between and is much larger than the above tolerance, we can conclude that -. Since , it appears that is a well-defined matrix of the full -rank.

Similarly, the matrix has the following non-zero singular values: while and . Hence, also is a well-defined matrix of the full -rank.

Let us represent in the form , where and are integers. In this case , i.e. , . Then, by the above-mentioned theorem 3.1 from [23], the left (right) partial indices of are stable if and only if , that is: (). Thus, we can conclude that -indices are and . Finally, since the matrices , are well-conditioned, we can expect that the partial indices of are also stable and equal to those computed numerically. But this needs anyway a rigorous proof.

Conclusion

We have used the explicit algorithm for the exact Wiener–Hopf factorization of a matrix polynomial with rational coefficients and with exact factorizable and described the main steps in its realization within the ExactMPF package. Its listing is included into the electronic supplementary material.

A message delivered by this paper is twofold: (a) the existence of an explicit algorithm does not necessarily guarantee that the factorization can be effectively performed; and, on the other hand, (b) in contrast to the general opinion, it is indeed possible to find sometimes exact values of the partial indices and then to perform exact factorization, even though the indices are not stable.

One of the main results of this paper is the necessary and sufficient condition for existence of the exact factorization for matrix polynomials whose coefficients are matrices with entries from the Gaussian rational field . It turns out that a matrix polynomial admits the exact Wiener–Hopf factorization if and only if the polynomial is exactly factorable.

The second (and extremely useful) result of the paper is the development of the package ExactMPF in the Maple software that realizes factorization algorithm exactly for a matrix polynomial with the coefficients from the Gaussian rational field and the exact factorizable determinant.

The use of the ExactMPF package is manyfold. The authors confidently declare that They also believe that it may be useful with additional efforts

it can be used for symbolic or symbolic-numeral factorization of matrix polynomials.

to analyse stability of the partial indices of matrix polynomials.

to make tests and experiments for verification of proved statements and for a motivation to formulate hypothesis or conjectures for further theoretical development of the theory.

in some cases to build an approximate factorization of a given matrix polynomial with arbitrary complex coefficients and

to solve real life problems as a computational tool. We shall demonstrate this in the accompanying paper dealing with solving of a discrete analogue of nonlinear Schrödinger equation by the inverse scattering transform method.