The Generalized Schur Algorithm and Some Applications

Abstract

The generalized Schur algorithm is a powerful tool allowing to compute classical decompositions of matrices, such as the $QR$ and $LU$ factorizations. When applied to matrices with particular structures, the generalized Schur algorithm computes these factorizations with a complexity of one order of magnitude less than that of classical algorithms based on Householder or elementary transformations. In this manuscript, we describe the main features of the generalized Schur algorithm. We show that it helps to prove some theoretical properties of the R factor of the $QR$ factorization of some structured matrices, such as symmetric positive definite Toeplitz and Sylvester matrices, that can hardly be proven using classical linear algebra tools. Moreover, we propose a fast implementation of the generalized Schur algorithm for computing the rank of Sylvester matrices, arising in a number of applications. Finally, we propose a generalized Schur based algorithm for computing the null-space of polynomial matrices.

1. Introduction

The generalized Schur algorithm (GSA) allows computing well-known matrix decompositions, such as $QR$ and $LU$ factorizations [1]. In particular, if the involved matrix is structured, i.e., Toeplitz, block-Toeplitz or Sylvester, the GSA computes the R factor of the $QR$ factorization with complexity of one order of magnitude less than that of the classical $QR$ algorithm [2], since it relies only on the knowledge of the so-called generators [2] associated to the given matrix, rather than on the knowledge of the matrix itself. The stability properties of the GSA are described in [3,4,5], where it is proven that the algorithm is weakly stable provided the involved hyperbolic rotations are performed in a stable way.

In this manuscript, we first show that, besides the efficiency properties, the GSA provides new theoretical insights on the bounds of the entries of the R factor of the $QR$ factorization of some structured matrices. In particular, if the involved matrix is a symmetric positive definite (SPD) Toeplitz or a Sylvester matrix, we prove that all or some of the diagonal entries of R monotonically decrease in absolute value.

We then propose a faster implementation of the algorithm described in [6] for computing the rank of a Sylvester matrix $S\in {\mathbb{R}}^{(m+n)\times (m+n)},$ whose entries are the coefficients of two polynomials of degree m and n, respectively. This new algorithm is based on the GSA for computing the R factor of the $QR$ factorization of $S.$ The proposed modification of the GSA-based method has a computational cost of $O(rl)$ floating point operations, where $l=min\{n,m\}$ and r is the computed numerical rank.

It is well known that the upper triangular factor R factor of the $QR$ factorization of a matrix $A\in {\mathbb{R}}^{n\times n}$ is equal to the upper triangular Cholesky factor $R_c\in {\mathbb{R}}^{n\times n}$ of $A^{T}A,$ up to a diagonal sign matrix D, i.e., $R=D{R}_c,$ $D=\mathrm{diag}(\pm 1,\cdots ,\pm 1)\in {\mathbb{R}}^{n\times n}.$ In this manuscript, we assume, without loss of generality, that the diagonal entries of R and $R_c$ are positive and since the matrices are then equal, we denote both matrices by R.

Finally, we propose a GSA-based approach for computing a null-space basis of a polynomial matrix, which is an important problem in several systems and control applications [7,8]. For instance, the computation of the null-space of a polynomial matrix arises when solving the column reduction problem of a polynomial matrix [9,10].

The manuscript is structured as follows. The main features of the GSA are provided in Section 2. In Section 3, a GSA implementation for computing the Cholesky factor R of a SPD Toeplitz matrix is described, which allows proving that the diagonal entries of R monotonically decrease. In Section 4, a GSA-based algorithm for computing the rank of a Sylvester matrix S is introduced, based on the computation of the Cholesky factor R of $S^{T}S.$ In addition, in this case, it is proven that the first diagonal entries of R monotonically decrease. The GSA-based method to compute the null-space of polynomial matrices is proposed in Section 5. The numerical examples are reported in Section 6 followed by the conclusions in Section 7.

2. The Generalized Schur Algorithm

Many of the classical factorizations of a symmetric matrix, e.g., $QR$ and $LD{L}^T,$ can be obtained by the GSA. If the matrix is Toeplitz-like, the GSA computes these factorizations in a fast way. For the sake of completeness, the basic concepts of the GSA for computing the R factor of the $QR$ factorization of structured matrices, such as Toeplitz and block-Toeplitz matrices, are introduced in this Section. A comprehensive treatment of the topic can be found in [1,2].

Let $A\in {\mathbb{R}}^{n\times n}$ be a symmetric positive definite (SPD) matrix. The semidefinite case is considered in Section 4 and Section 5. The displacement of A with respect to a matrix Z of order $n,$ is defined as

$${\nabla }_{Z}A=A-ZA{Z}^T,$$

(1)

while the displacement rank k of A with respect to Z is defined as the rank of ${\nabla }_{Z}A.$ If $\mathrm{rank}({\nabla }_{Z}A)=k,$ Equation (1) can be written as the sum of k rank-one matrices,

$${\nabla }_{Z}A=\sum _{i=1}^{k_1}{\mathbf{g}}_i^{(p)}{{\mathbf{g}}_i^{(p)}}^T-\sum _{i=1}^{k_2}{\mathbf{g}}_i^{(n)}{{\mathbf{g}}_i^{(n)}}^T,$$

where $(k_1,n-k_1-k_2,k_2)$ is the inertia of ${\nabla }_{Z}A,$ $k=k_1+k_2,$ and the vectors ${\mathbf{g}}_i^{(p)}\in {\mathbb{R}}^n,i=1,\dots ,k_1,$ ${\mathbf{g}}_i^{(n)}\in {\mathbb{R}}^n,i=1,\dots ,k_2,$ are called the positive and the negative generators of A with respect to $Z,$ respectively, conversely, if there is no ambiguity, simply the positive and negative generators of A. The matrix $G\equiv {[{\mathbf{g}}_1^{(p)},{\mathbf{g}}_2^{(p)},\dots ,{\mathbf{g}}_{k_1}^{(p)},{\mathbf{g}}_1^{(n)},{\mathbf{g}}_2^{(n)},\dots ,{\mathbf{g}}_{k_2}^{(n)}]}^T$ is called the generator matrix.

The matrix Z is a nilpotent matrix. In particular, for Toeplitz and block-Toeplitz matrices, the matrix Z can be chosen as the shift and the block shift matrix

$$Z_1=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ 1& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & 1& 0\end{array}\right],Z_2=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ Z_1& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & Z_1& 0\end{array}\right],$$

respectively.

The implementation of the GSA relies only on the knowledge of the generators of A rather than on the knowledge of the matrix itself [1].

Let

$$J=\mathrm{diag}(\underset{k_1}{\underbrace{1,1,\dots ,1}},\underset{k_2}{\underbrace{-1,-1,\dots ,-1}}).$$

Since

$$\begin{array}{ccc}A-ZA{Z}^T\hfill & =& G^{T}JG,\hfill \\ ZA{Z}^T-Z^{2}A{Z^2}^T\hfill & =& Z{G}^{T}JG{Z}^T,\hfill \\ ⋮\hfill & & ⋮\hfill \\ Z^{n-2}A{Z^{n-2}}^T-Z^{n-1}A{Z^{n-1}}^T\hfill & =& Z^{n-2}{G}^{T}JG{Z^{n-2}}^T,\hfill \\ Z^{n-1}A{Z^{n-1}}^T\hfill & =& Z^{n-1}{G}^{T}JG{Z^{n-1}}^T,\hfill \end{array}$$

(2)

then, adding all members of the left and right-hand sides of Equation (2) yields

$$A=\sum _{j=0}^{n-1}{Z}^j{G}^{T}JG{Z^j}^T,$$

(3)

which expresses the matrix A in terms of its generators.

Exploiting Equation (2), we show how the GSA computes R by describing its first iteration. Observe that the matrix products involved in the right-hand side of Equation (2) have their first row equal to zero, with the exception of the first product, $G^{T}JG$.

A key role in GSA is played by J-orthogonal matrices [11,12], i.e., matrices $\Phi$ satisfying ${\Phi }^{T}J\Phi =J.$

Any such matrix $\Phi$ can be constructed in different ways [11,12,13,14]. For instance, it can be considered as the product of Givens and hyperbolic rotations. In particular, a Givens rotation acting on rows i and j of the generator matrix is chosen if $J(i,i)J(j,j)>0,$ $i,j\in \{1,\dots ,n\},i\ne j.$ Otherwise, a hyperbolic rotation is considered. Indeed, suitable choices of $\Phi$ allow efficient implementations of GSA, as shown in Section 4.

Let $G_0\equiv G$ and ${\Phi }_1$ be a J-orthogonal matrix such that

$${\tilde{G}}_1={\Phi }_1{G}_0,{\tilde{G}}_1{\mathbf{e}}_1={[{\alpha }_1,0,\dots ,0]}^T,\mathrm{with}{\alpha }_1>0,$$

(4)

and ${\mathbf{e}}_i,i=1,\dots ,n,$ be the ith column of the identity matrix. Furthermore, let ${\tilde{\mathbf{g}}}_1^T$ and ${\tilde{\Gamma }}_1$ be the first and last $k-1$ rows of ${\tilde{G}}_1,$ respectively, i.e., ${\tilde{G}}_1=\left[\begin{array}{c}{\tilde{\mathbf{g}}}_1^T\\ {\tilde{\Gamma }}_1\end{array}\right].$

From Equation (4), it turns out that the first column of ${\tilde{\Gamma }}_1$ is zero. Let $\tilde{J}$ be the matrix obtained by deleting the first row and column from J. Then, Equation (2) can be written as follows,

$$\begin{array}{ccc}\hfill A& =& \sum _{j=0}^{n-1}{Z}^j{G}_0^{T}J{G}_0{Z^j}^T\hfill \\ \hfill & =& \sum _{j=0}^{n-1}{Z}^j{G}_0^T{\Phi }_1^{T}J{\Phi }_1{G}_0{Z^j}^T\hfill \\ \hfill & =& \sum _{j=0}^{n-1}{Z}^j{\left[\begin{array}{c}{\tilde{\mathbf{g}}}_1^T\\ {\tilde{\Gamma }}_1\end{array}\right]}^{T}J\left[\begin{array}{c}{\tilde{\mathbf{g}}}_1^T\\ {\tilde{\Gamma }}_1\end{array}\right]{Z^j}^T\hfill \\ \hfill & =& {\tilde{g}}_1{\tilde{g}}_1^T+\sum _{j=1}^{n-1}{Z}^j{\tilde{g}}_1{\tilde{g}}_1^T{Z^j}^T+\sum _{j=0}^{n-2}{Z}^j{\tilde{\Gamma }}_1^T\tilde{J}{\tilde{\Gamma }}_1{Z^j}^T+\underset{=0}{\underbrace{Z^{n-1}{\tilde{\Gamma }}_1^T\tilde{J}{\tilde{\Gamma }}_1{Z^{n-1}}^T}}\hfill \\ \hfill & =& {\tilde{\mathbf{g}}}_1{\tilde{\mathbf{g}}}_1^T+\sum _{j=0}^{n-2}{Z}^j{\left[\begin{array}{c}{\tilde{\mathbf{g}}}_1^T{Z}^T\\ {\tilde{\Gamma }}_1\end{array}\right]}^{T}J\left[\begin{array}{c}{\tilde{\mathbf{g}}}_1^T{Z}^T\\ {\tilde{\Gamma }}_1\end{array}\right]{Z^j}^T\hfill \\ \hfill & =& {\tilde{\mathbf{g}}}_1{\tilde{\mathbf{g}}}_1^T+\sum _{j=0}^{n-2}{Z}^j{G}_1^{T}J{G}_1{Z^j}^T,\hfill \\ \hfill & =& {\tilde{\mathbf{g}}}_1{\tilde{\mathbf{g}}}_1^T+A_1,\hfill \end{array}$$

where $G_1\equiv {[Z{\tilde{\mathbf{g}}}_1,{\tilde{\Gamma }}_1^T]}^T,$ that is, $G_1$ is obtained from ${\tilde{G}}_1$ by multiplying ${\tilde{g}}_1$ with $Z,$ and $A_1\equiv {\sum }_{j=0}^{n-2}{Z}^j{G}_1^{T}J{G}_1{Z^j}^T.$ If A is a Toeplitz matrix, this multiplication with Z corresponds to displacing the entries of ${\tilde{g}}_1$ one position downward, while it corresponds to a block-displacement downward in the first generator if A is a block-Toeplitz matrix.

Thus, the first column of $G_1$ is zero and, hence, ${\tilde{\mathbf{g}}}_1^T$ is the first row of the R factor of the $QR$ factorization of $A.$ The above procedure is recursively applied to $A_1$ to compute the other rows of $R.$

The jth iteration of GSA, $j=1,\dots ,n,$ involves the products ${\Phi }_j{G}_{j-1}$ and $Z{\tilde{\mathbf{g}}}_1$. The former multiplication can be computed in $O\left(k(n-j)\right)$ operations [11,12], and the latter is done for free if Z is either a shift or a block–shift matrix. Therefore, if the displacement rank k of A is small compared to $n,$ the GSA computes the R factor in $O(k{n}^2)$ rather than in $O(n^3)$ operations, as required by standard algorithms [15].

For the sake of completeness, the described GSA implementation is reported in the following matlab style function. (The function givens is the matlab function having as input two scalars, $x_1$ and $x_2,$ and as output an orthogonal $2\times 2$ matrix $\Theta$ such that $\Theta \left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}\sqrt{x_1^2+x_2^2}\\ 0\end{array}\right].$ The function Hrotate computes the coefficients of the $2\times 2$ hyperbolic rotation $\Phi$ such that, given two scalars $x_1$ and $x_2,$ $|x_1|>|x_2|,$ $\Phi \left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}\sqrt{x_1^2-x_2^2}\\ 0\end{array}\right].$ The function Happly applies $\Phi$ to two rows of the generator matrix. Both functions are defined in [12]).

function$[R]=$GSA$(G,n);$

fori = 1 : n,

  for j = 2 : k₁,

    $\Theta =$ givens $(G(1,i),G(j,i));$

    $G([1,j],i:n)=\Theta \ast G([1,j],i:n);$

  end % for

  for j = k₁ + 2 : k₁ + k₂,

    Θ = govens(G(k₁ + 1, i), G(j, i));

    G([k₁ + 1, j], i : n) = Θ ∗ G[k₁ + 1, j], i : n);

  end % for

  [c₁, s₁] = Hrotate(G(1, i), G(k₁ + 1, i));

  G([1, k₁ + 1], i : n) = Happly(c₁, s₁, G([1, k₁ + 1], i : n), n − i + 1);

  R(i, i : n) = G(1, i : n);

  G(1, i + 1 : n) = G(1, i : n − 1); G(1, i) = 0;

end % for

The GSA has been proven to be weakly stable [3,4], provided the hyperbolic transformations involved in the construction of the matrices ${\Phi }_j$ are performed in a stable way [3,11,12].

3. GSA for SPD Toeplitz Matrices

In this section, we describe the GSA for computing the R factor of the Cholesky factorization of a SPD Toeplitz matrix A, with R upper triangular, i.e., $A=R^{T}R.$ Moreover, we show that the diagonal entries of R decrease monotonically.

Let $A\in {\mathbb{R}}^{n\times n}$ and $Z\in {\mathbb{R}}^{n\times n}$ be a SPD Toeplitz matrix and a shift matrix, respectively, i.e.,

$$A=\left[\begin{array}{cccc}{t}_1& t_2& \ddots & t_n\\ t_2& \ddots & \ddots & \ddots \\ \ddots & \ddots & \ddots & t_2\\ t_n& \ddots & t_2& t_1\end{array}\right],Z_n=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ 1& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & 1& 0\end{array}\right],$$

and let $\mathbf{t}=A(:,1).$ Then,

$${\nabla }_{Z}A=\left[\begin{array}{cccc}{t}_1& t_2& \cdots & t_n\\ t_2& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮\\ t_n& 0& \cdots & 0\end{array}\right],$$

i.e., ${\nabla }_{Z}A$ is a symmetric rank-2 matrix. Moreover, the generator matrix G is given by

$$G=\left[\begin{array}{c}{\mathbf{g}}_1^T\\ {\mathbf{g}}_2^T\end{array}\right],\mathrm{with}{\mathbf{g}}_1=\frac{\mathbf{t}}{\sqrt{t_1}},{\mathbf{g}}_2={[0,{\mathbf{g}}_1{(2:n)}^T]}^T.$$

In this case, the GSA can be implemented in matlab-like style as follows.

function$[R]=$GSA_chol$(G_0)$

fori = 1 : n,

  [c₁, s₁] = Hrotate(G_{i − 1}(1, i), G⁽ⁱ⁾(2, i)); G_{i − 1}(:, i : n) = Happly(c₁, s₁, G_{i − 1}(:, i : n), n − i + 1);

  R(i, i : n) = G_{i − 1}(1, i : n);

  G_i(1, i + 1 : n) = G_{i − 1}(1, i : n − 1); G_i(2, i + 1 : n) = G_{i − 1}(2, i + 1 : n − 1);

end % for

The following lemma holds.

Lemma 1.

Let A be a SPD Toeplitz matrix and let R be its Cholesky factor, with R upper triangular. Then,

$$R(i-1,i-1)\ge R(i,i),i=2,\dots ,n.$$

Proof. 

At each step i of GSA_chol, $i=1,\dots ,n,$ first a hyperbolic rotation is applied to $G_{i-1}$ in order to annihilate the element $G_i(2,i).$ Hence, the first row of $G_{i-1}$ becomes the row i of $R.$ Finally, $G_i(1,:)$ is obtained displacing the entries of the first row of $G_{i-1}$ one position right, while $G_i(2,:)$ is equal to $G_{i-1}(2,:).$ Taking into account that $G_{i-1}(2,1)=0,$ the diagonal entries of R are

$$\begin{array}{ccc}\hfill R(1,1)& =& G_0(1,1)\hfill \\ \hfill R(2,2)& =& \sqrt{G_1^2(1,2)-G_1^2(2,2)}=\sqrt{R^2(1,1)-G_1^2(2,2)}\le R(1,1);\hfill \\ \hfill & ⋮& \hfill \\ \hfill R(i,i)& =& \sqrt{G_{i-1}^2(1,i)-G_{i-1}^2(2,i)}=\sqrt{R^2(i-1,i-1)-G_{i-1}^2(2,i)}\le R(i-1,i-1);\hfill \\ \hfill & ⋮& \hfill \\ \hfill R(n,n)& =& \sqrt{G_{n-1}^2(1,n)-G_{n-1}^2(2,n)}=\sqrt{R^2(n-1,n-1)-G_{n-1}^2(2,n)}\le R(n-1,n-1).\hfill \end{array}$$

 □

4. Computing the Rank of Sylvester Matrices

In this section, we focus on the computation of the rank of Sylvester matrices. The numerical rank of a Sylvester matrix is a useful information for determining the degree of the greatest common divisor of the involved polynomials [6,16,17].

A GSA-based algorithm for computing the rank of S has been recently proposed in [6]. It is based on the computation of the Cholesky factor R of $S^{T}S,$ with R upper triangular, i.e., $R^{T}R=S^{T}S.$

Here, we propose a more efficient variant of this algorithm that allows proving that the first entries of R monotonically decrease.

Let $w_i\in \mathbb{R},i=0,1,\dots ,n,$ and let $y_i\in \mathbb{R},i=0,1,\dots ,m.$ Denote by $w(x)$ and $y(x)$ two univariate polynomials,

$$\begin{array}{cc}w(x)=w_n{x}^n+w_{n-1}{x}^{n-1}+\cdots +w_{1}x+w_0,\hfill & w_n\ne 0,\hfill \\ y(x)=y_m{x}^m+y_{m-1}{x}^{m-1}+\cdots +y_{1}x+y_0,\hfill & y_m\ne 0.\hfill \end{array}$$

(5)

Let $S\in {\mathbb{R}}^{(m+n)\times (m+n)}$ be the Sylvester matrix defined as follows,

$$S=\left[\begin{array}{cc}W& Y\end{array}\right],W=\left[\begin{array}{cccc}{w}_n& & & \\ w_{n-1}& w_n& & \\ ⋮& w_{n-1}& \ddots & \\ w_1& ⋮& \ddots & w_n\\ w_0& w_1& \ddots & w_{n-1}\\ & w_0& \ddots & ⋮\\ & & \ddots & w_1\\ & & & w_0\end{array}\right],Y=\left[\begin{array}{cccc}{y}_m& & & \\ y_{m-1}& y_m& & \\ ⋮& y_{m-1}& \ddots & \\ y_1& ⋮& \ddots & y_m\\ y_0& y_1& \ddots & y_{m-1}\\ & y_0& \ddots & ⋮\\ & & \ddots & y_1\\ & & & y_0\end{array}\right],$$

(6)

with $W\in {\mathbb{R}}^{(m+n)\times m}$ and $Y\in {\mathbb{R}}^{(m+n)\times n}$ band Toeplitz matrices.

We now describe how the GSA-based algorithm proposed in [6] for computing the rank of S can be implemented in a faster way. This variant is based on the computation of the Cholesky factor $R\in {\mathbb{R}}^{(m+n)\times (m+n)}$ of $S^{T}S,$ with R upper triangular, i.e., $R^{T}R=S^{T}S.$

Defining

$$Z=\left[\begin{array}{c}{Z}_m\\ & Z_n\end{array}\right],\mathrm{with}{Z}_k={\left[\begin{array}{cccc}0& 0& \cdots & 0\\ 1& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & ⋮\\ 0& \cdots & 1& 0\end{array}\right]}_{k\times k},k\in \mathbb{N},$$

(7)

the generator matrix G of $S^{T}S$ with respect to Z is then given by [6]

$$G={\left[\begin{array}{cccc}{\mathbf{g}}_1& {\mathbf{g}}_2& {\mathbf{g}}_3& {\mathbf{g}}_4\end{array}\right]}^T$$

where

$$\begin{array}{c}{\mathbf{g}}_1={\mathbf{x}}_1/{\parallel S(:,1)\parallel }_2,\hfill \\ {\mathbf{g}}_2([2:n+m])={\mathbf{x}}_2([2:n+m])/{\parallel S(:,m+1)\parallel }_2,{\mathbf{g}}_2(1)=0,\hfill \\ {\mathbf{g}}_3(2:n+m)={\mathbf{g}}_1(2:n+m),{\mathbf{g}}_3(1)=0,\hfill \\ {\mathbf{g}}_4([1:m,m+2:n+m])={\mathbf{g}}_2([1:m,m+2:n+m]),{\mathbf{g}}_4(m+1)=0,\hfill \end{array}$$

(8)

with ${\mathbf{x}}_1=S^{T}S{\mathbf{e}}_1,$ ${\mathbf{x}}_2=S^{T}S{\mathbf{e}}_{m+1},$ ${\mathbf{e}}_j$ the jth vector of the canonical basis of ${\mathbb{R}}^{m+n},$ and $J=\mathrm{diag}(1,1,-1,-1)$.

The algorithm proposed in [6] is based on the following GSA implementation for computing the R factor of the $QR$ factorization of $S.$

function$[R]=$GSA_chol2$(G)$

fori = 1 : n,

  Θ₁ = givens(G(1, i), G(2, i));    Θ₂ = givens(G(3, i), G(4, i));

  G(1 : 2, i : n) = Θ₁G(1 : 2, i : n);   G(3 : 4, i : n) = Θ₂G(3 : 4, i : n);

  [c₁, s₁] = Hrotate(G(1, i), G(3, i));

  G([1,3], i : n) = Happly(c₁, s₁, G([1,3], i : n), n − i + 1);

  R(i, i : n) = G_i(1, i : n);

  G(1, i + 1 : n) = G(1, i : n − 1)Z^T;   G(2, i + 1 : n) = G(2 : 4, i + 1 : n − 1);

end % for

At the ith iteration of the algorithm, $i=1,\dots ,n,$ the G_ivens rotations ${\Theta }_1$ and ${\Theta }_2$ are computed and applied, respectively, to the first and second generators, and to the third and fourth generators,   to annihilate $G(2,i)$ and $G(4,i).$ Hence, the hyperbolic rotation $\left[\begin{array}{cc}{c}_1& -s_1\\ -s_1& c_1\end{array}\right]$ is applied to the first and the third row of G to annihilate $G(3,i).$ Finally, the first row of G becomes the ith row of R and the first row of G is multiplied by $Z^T.$

Summarizing, at the first step of the ith iteration of GSA, all entries of the ith column but the first one of G, are annihilated. If the number of rows of G is greater than $2,$ this can be accomplished in different ways (see [5,14]).

Analyzing the pattern of the generators in Equation (8), we are able to derive a different implementation of $GSA$ that costs $O(rl)$, with $l=min\{n,m\}$. Moreover, this implementation allows proving that the first l diagonal entries of R are monotonically decreasing.

We observe that the matrix $W^{T}W$ in Equation (6) is the SPD Toeplitz matrix

$$W^{T}W={\left[\begin{array}{ccccccc}{t}_1& t_2& \cdots & t_n& t_{n+1}& & \\ t_2& t_1& t_2& \ddots & t_n& \ddots & \\ ⋮& t_2& \ddots & \ddots & \ddots & \ddots & t_{n+1}\\ t_n& ⋮& \ddots & \ddots & \ddots & ⋮& t_n\\ t_{n+1}& t_n& \ddots & \ddots & \ddots & t_2& ⋮\\ & \ddots & \ddots & ⋮& t_2& t_1& t_2\\ & & t_{n+1}& t_n& \cdots & t_2& t_1\end{array}\right]}_{m\times m},$$

(9)

with

$$t_i=\sum _{j=i}^{n+1}{w}_{j-1}{w}_{j-i},i=1,2,\dots ,n+1.$$

Since

$$S^{T}S=\left[\begin{array}{cc}{W}^{T}W& W^{T}Y\\ Y^{T}W& Y^{T}Y\end{array}\right],$$

if $n\ll m,$ from Equation (9), it turns out that $G([1,3],n+2:m)=0.$ Moreover, the rows $G(2,:)$ and $G(4,:)$ have their first entry equal to zero and differ only in their entry in column $m+1.$ This particular pattern of G is close to the ones described in [13,14,18], allowing to design an alternative GSA implementation with respect to that considered in [6], and thereby reducing the complexity from $O(r(n+m))$ to $O(rl)$, where r is the computed rank of S and $l=min\{n,m\}.$

Since the description of the above GSA implementation is quite cumbersome and similar to the algorithms reported in [13,14,18], we omit it here. The corresponding matlab pseudo–code can be obtained from the authors upon request.

If the matrix S has rank $r<(n+m),$ at the $k=(n+m-r+1)$st iteration, it turns out that $G^2(1,k)-G^2(3,k)=0$ in exact arithmetic [6]. Therefore, at each iteration of the algorithm we check whether

$$G^2(1,k)-G^2(3,k)>tol,$$

(10)

where $tol$ is a fixed tolerance. If Equation (10) is not satisfied, we stop the computation considering k as the computed numerical rank of $S.$

The R factor of the $QR$ factorization of S is unique if the diagonal entries of R are positive. The considered GSA implementation, yielding the rank of S and based on computing the R factor of the $QR$ factorization of S, allows us to prove that the first l entries of the diagonal of R are ordered in a decreasing order, with $l=min\{m,n\}.$ In fact, the following theorem holds.

Theorem 1.

Let $R^{T}R=S^{T}S$ be the Cholesky factorization of $S^{T}S$ with S the Sylvester matrix defined in Equation (6) with rank $r\ge l=min\{m,n\}.$ Then,

$$R(i-1,i-1)\ge R(i,i)\ge 0,i=2,\dots ,l.$$

(11)

Proof. 

Each entry i of the diagonal of R is determined by the ith entry of the first row of G at the end of iteration $i,$ for $i=1,\dots ,m+n.$ Let us define $\widehat{G}\equiv G(:,1:l)$ and consider the following alternative implementation of the GSA for computing the first l columns of the Cholesky factor of $S^{T}S.$

for$i=1:l,$

  $\Theta =$ givens$(\widehat{G}(1,i),\widehat{G}(2,i));$

  $\widehat{G}(1:2,i:l)=\Theta \ast \widehat{G}(1:2,i:l);$

  $[c_1,s_1]=$ Hrotate$(\widehat{G}(1,i),\widehat{G}(4,i));$ $\widehat{G}([1,4],:)=\mathtt{Happly}(c_1,s_1,\widehat{G}([1,4],:),l);$

  $[c_2,s_2]=$ Hrotate$(\widehat{G}(1,i),\widehat{G}(3,i));$ $\widehat{G}([1,3],:)=\mathtt{Happly}(c_2,s_2,\widehat{G}([1,3],:),l);$

  $R(i,i:l)=\widehat{G}(1,i:l);$

  $\widehat{G}(1,i+1:l)=\widehat{G}(1,i:l-1);$

  $\widehat{G}(1,i)=0;$

end % for

We observe that, for $i=1$, $\widehat{G}(1,1)$ is the only entry in the first column of $\widehat{G}$ different from $0.$ Hence, $R(1,i)=\widehat{G}(1,1:l)$ and the first iteration amounts only to shifting $\widehat{G}(1,1:l)$ one position rightward, i.e., $\widehat{G}(1,2:l)=\widehat{G}(1,1:l-1),\widehat{G}(1,1)=0.$

At the beginning of iteration $i=2,$ the second and the fourth row of $\widehat{G}$ are equal Equation (8). Hence, when applying a G_ivens rotation to the first and the second row in order to annihilate the entry $\widehat{G}(2,i)$ and when subsequently applying a hyperbolic rotation to the first and fourth row of $\widehat{G}$ in order to annihilate $\widehat{G}(4,i),$ it turns out that $\widehat{G}(2,i:l)$ and $\widehat{G}(4,i:l)$ are then modified but still equal to each other, while $\widehat{G}(1,i:l)$ remains unchanged. The equality between $\widehat{G}(2,:)$ and $\widehat{G}(4,:)$ is maintained throughout the iterations $1,2,\dots ,l.$

Therefore, the second and the fourth row of $\widehat{G}$ do not play any role in computing $R(1:l,1:l)$ and can be neglected. Hence, the GSA for computing $R(1:l,1:l)$ reduces only to applying a hyperbolic rotation to the first and the third generators, as described in the following algorithm.

for$i=1:l,$

  $[c_2,s_2]=$ Hrotate$(\widehat{G}(1,i),\widehat{G}(3,i));$ $\widehat{G}([1,3],:)=\mathtt{Happly}(c_2,s_2,\widehat{G}([1,3],:),l);$

  $R(i,i:l)=\widehat{G}(1,i:l);$

  $\widehat{G}(1,i+1:l)=\widehat{G}(1,i:l-1);$

  $\widehat{G}(1,i)=0;$

end % for

Since at the beginning of iteration $i,i=2,\dots ,i,$ $\widehat{G}(1,i:l)=R(i-1,i-1:l-1),$ then the involved hyperbolic rotation $\Phi =\left[\begin{array}{cc}{c}_2& -s_2\\ -s_2& c_2\end{array}\right]$ is such that

$$\Phi \left[\begin{array}{c}\widehat{G}(1,i)\\ \widehat{G}(3,i)\end{array}\right]=\Phi \left[\begin{array}{c}R(i-1,i-1)\\ \widehat{G}(3,i)\end{array}\right]=\left[\begin{array}{c}\widehat{G}(1,i)\\ 0\end{array}\right]=\left[\begin{array}{c}R(i,i)\\ 0\end{array}\right],$$

where the updated $\widehat{G}(1,i)$ is equal to $\sqrt{\widehat{G}{(1,i)}^2-\widehat{G}{(3,i)}^2}\ge 0.$ Therefore, $R(i,i)=\sqrt{R{(i-1,i-1)}^2-\widehat{G}{(3,i)}^2}\ge 0,$ and thus $R(i,i)\le R(i-1,i-1).$ □

Remark 1.

The above GSA implementation allows to prove the inequality Equation (11). This property is difficult to obtain if the $QR$ factorization is performed via Householder transformations or if the classical Cholesky factorization of $S^{T}S$ is used.

5. GSA for Computing the Null-Space of Polynomial Matrices

In this section, we consider the problem of computing a polynomial basis $X(s)\in {\mathbb{R}}^{n\times (n-\rho )}$ of the null-space of an $m\times n$ polynomial matrix of degree $\delta$ and rank $\rho \le min(m,n),$

$$M(s)=\sum _{i=0}^{\delta }{M}_i{s}^i,M_i\in {\mathbb{R}}^{m\times n},i=0,\dots ,\delta .$$

(12)

As described in [8,19,20], the above problem is equivalent to that of computing the null-space of a related block-Toeplitz matrix. Algorithms to solve this problem are proposed in [8,19] but they do not explicitly exploit the structure of the involved matrix. Algorithms to solve related problems have also been described in the literature, e.g., in [8,19,21,22].

In this paper, we propose an algorithm for computing the null-space of polynomial matrices based on a variant of the GSA for computing the null-space of a related band block-Toeplitz matrix [8].

5.1. Null-Space of Polynomial Matrices

A polynomial vector $\mathbf{v}(s)={\sum }_{i=0}^{\gamma }{\mathbf{v}}_i{s}^i,{\mathbf{v}}_i\in {\mathbb{R}}^n,i=0,\dots ,\gamma ,$ $\gamma \in \mathbb{N},$ is said to belong to the null-space of (12) if

$$M(s)\mathbf{v}(s)=0\iff \sum _{j=0}^{\delta }{M}_j{s}^j\sum _{i=0}^{\gamma }{\mathbf{v}}_i{s}^i=0.$$

The polynomial vector $\mathbf{v}(s)$ belongs to the null-space of $M(s)$ iff $\mathbf{v}={[{\mathbf{v}}_0^T,{\mathbf{v}}_1^T,\dots ,{\mathbf{v}}_{\gamma }^T]}^T,$ ${\mathbf{v}}_i\in {\mathbb{R}}^n,i=0,\dots ,\gamma ,$ is a vector belonging to the null-space of the band block-Toeplitz matrix

$$T={\left[\begin{array}{ccccc}{M}_0& & & & \\ M_1& M_0& & & \\ ⋮& M_1& \ddots & & \\ M_{\delta }& ⋮& \ddots & M_0& \\ & M_{\delta }& \ddots & M_1& M_0\\ & & \ddots & ⋮& M_1\\ & & & M_{\delta }& ⋮\\ & & & & M_{\delta }\end{array}\right]}_{\widehat{m}\times \widehat{n}},$$

(13)

where $\widehat{m}=m(\delta +n_b),\widehat{n}=n{n}_b,$ with $n_b=\gamma +1$ the number of block columns of $T,$ that can be determined, e.g., by the algorithm described in [8]. Hence, the problem of computing the null-space of the polynomial matrix in Equation (12) is equivalent to the problem of computing the null-space of the matrix in Equation (13). To normalize the entries in this matrix, it is appropriate to first perform a QR factorization of each block column of T:

$$\left[\begin{array}{c}{M}_0\\ M_1\\ ⋮\\ M_{\delta }\end{array}\right]=\left[\begin{array}{c}{Q}_0\\ Q_1\\ ⋮\\ Q_{\delta }\end{array}\right]U,\mathrm{where}\sum _i{Q}_i^T{Q}_i=I_n,$$

and to absorb the upper triangular factor U in the vector $\mathbf{u}(s):=U\mathbf{v}(s)$. The convolution equation $M(s)\mathbf{v}(s)=0$ then becomes an equation of the type $Q(s)\mathbf{u}(s)=0,$ but where the coefficient matrices $Q_i$ of $Q(s)$ form together an orthonormalized matrix.

Remark 2.

Above, we have assumed that there are no constant vectors $\mathbf{v}$ in the kernel of $M(s)$. If there are, then, the block column of $M_i$ matrices has rank less than n and the above factorization will discover it in the sense that the matrix U is nonsquare and the matrices $Q_i$ have less columns than $M_i,i=0,1,\dots ,\delta .$ This trivial null-space can be eliminated and we therefore assume that the rank was full. For simplicity, from now on, we also assume that the coefficient matrices of the polynomial matrix $M(s)$ were already normalized in this way and the norm of the block columns of T are thus orthonormalized. This normalization proves to be very useful in the sequel.

Denote by

$$Z_{n_b}={\left[\begin{array}{cccc}{0}_n& & & \\ I_n& 0_n& & \\ & \ddots & \ddots & \\ & & I_n& 0_n\end{array}\right]}_{\widehat{n}\times \widehat{n}}\mathrm{and}Z={\left[\begin{array}{cc}{Z}_{n_b}& \\ & Z_{n_b}\end{array}\right]}_{2\widehat{n}\times 2\widehat{n}},$$

where $0_n$ is the null–matrix of order $n\in \mathbb{N}.$

If ${\mathbf{v}}_i\ne 0,$ $0<i<\gamma ,$ and ${\mathbf{v}}_j=0,j=i+1,\dots ,\gamma ,$ i.e.,

$$\mathbf{v}={[{\mathbf{v}}_0^T,{\mathbf{v}}_1^T,\dots ,{\mathbf{v}}_i^T,\underset{\gamma -i}{\underbrace{0,\dots ,0}}]}^T,$$

and $\mathbf{v}\in ker(T),$ then also $Z_{n_b}^k\mathbf{v}\in ker(T),k=0,1,\dots ,\gamma -i.$ In this case, the vector $\mathbf{v}$ is said to be a generator vector of a chain of length $\gamma -i+1$ of the null-space of $T.$

The proposed algorithm for the computation of the null-space of polynomial matrices is based on the GSA for computing the R factor of the $QR$-factorization of the matrix T in Equation (13) and, if R is full column rank, its inverse $R^{-1}$.

Let us first assume that the matrix T is full rank, i.e., $\mathrm{rank}(T)=\rho =min\{\widehat{m},\widehat{n}\}.$ Without loss of generality, we suppose $\widehat{m}\ge \widehat{n}.$ If $\widehat{m}<\widehat{n},$ the algorithm still computes the R factor in trapezoidal form [23]. Moreover, in this case, we compute the first $\widehat{m}$ rows of the inverse of the matrix obtained appending the last $\widehat{n}-\widehat{m}$ rows of the identity matrix of order $\widehat{n}$ to $R.$

Let us consider the SPD block-Toeplitz matrix

$$\widehat{T}=T^{T}T=\left[\begin{array}{cccccc}{\widehat{T}}_0& {\widehat{T}}_1& \cdots & {\widehat{T}}_{\delta }& & \\ {\widehat{T}}_1& {\widehat{T}}_0& {\widehat{T}}_1& \ddots & \ddots & \\ ⋮& {\widehat{T}}_1& \ddots & \ddots & \ddots & {\widehat{T}}_{\delta }\\ {\widehat{T}}_{\delta }& \ddots & \ddots & \ddots & \ddots & ⋮\\ & \ddots & \ddots & \ddots & {\widehat{T}}_0& {\widehat{T}}_1\\ & & {\widehat{T}}_{\delta }& \cdots & {\widehat{T}}_1& {\widehat{T}}_0\end{array}\right]\in {\mathbb{R}}^{\widehat{n}\times \widehat{n}},$$

(14)

whose blocks are

$${\widehat{T}}_{i-j}=\left\{\begin{array}{cc}\sum _{k=0}^{\delta -|i-j|}{M}_{k+|i-j|}^T{M}_k,\hfill & \mathrm{if}i\le j\hfill \\ \sum _{k=0}^{\delta -|i-j|}{M}_k^T{M}_{k+|i-j|},\hfill & \mathrm{if}i>j.\hfill \end{array}\right.$$

(15)

Notice that, because of the normalization introduced before, we have that ${\widehat{T}}_0=I_n$ and $\parallel {\widehat{T}}_i{\parallel }_2\le 1$. This is used below. The matrix

$$W=\left[\begin{array}{cc}\widehat{T}& I_{\widehat{n}}\\ I_{\widehat{n}}& 0_{\widehat{n}}\end{array}\right]$$

(16)

can be factorized in the following way,

$$W={\widehat{R}}^T\widehat{J}\widehat{R}\equiv \left[\begin{array}{cc}{R}^T& \\ R^{-1}& R^{-1}\end{array}\right]\left[\begin{array}{cc}{I}_{\widehat{n}}& \\ & -I_{\widehat{n}}\end{array}\right]\left[\begin{array}{cc}R& R^{-T}\\ & R^{-T}\end{array}\right],$$

(17)

where $R\in {\mathbb{R}}^{\widehat{n}\times \widehat{n}}$ is the factor R of the $QR$-factorization of $T,$ i.e., the Cholesky factor of $\widehat{T}.$ Hence, R and its inverse $R^{-1}$ can be retrieved from the first $\widehat{n}$ columns of the matrix ${\widehat{R}}^T.$

The displacement matrix and the displacement rank of W with respect to $Z,$ are given by

$${\nabla }_Z(W)=W-ZW{Z}^T=\left[\begin{array}{ccccccccccc}\hfill I_n& \hfill {\widehat{T}}_1& \hfill \cdots & {\widehat{T}}_{\delta }& 0_n& \cdots & 0_n& I_n& 0_n& \cdots & 0_n\\ \hfill {\widehat{T}}_1\\ \hfill ⋮\\ \hfill {\widehat{T}}_{\delta }\\ \hfill 0_n\\ \hfill ⋮\\ \hfill 0_n\\ \hfill I_n\\ \hfill 0_n\\ \hfill ⋮\\ \hfill 0_n\end{array}\right]$$

(18)

and $\rho (W,Z)=\mathrm{rank}({\nabla }_Z(W)),$ respectively, with ${\nabla }_Z(W)\in {\mathbb{R}}^{2\widehat{n}\times 2\widehat{n}}.$

Then, taking the order n of the matrices ${\widehat{T}}_i,i=0,1,\dots ,\delta ,$ into account, it turns out that $\rho (W,Z)\le 2n.$

Hence, Equation (18) can be written as the difference of two matrices of rank at most n, i.e.,

$${\nabla }_Z(W)={G^{(+)}}^T{G}^{(+)}-{G^{(-)}}^T{G}^{(-)}=G^{T}JG,\mathrm{where}G:=\left[\begin{array}{c}{G}^{(+)}\\ G^{(-)}\end{array}\right]\mathrm{and}J=\mathrm{diag}(I_n,-I_n).$$

Since ${\widehat{T}}_0=I_n,$ the construction of G does not require any computation: it is easy to check that G is given by

$$G:=\left[\begin{array}{c}{G}^{(+)}\\ G^{(-)}\end{array}\right]=\left[\begin{array}{ccccccccccc}{I}_n& {\widehat{T}}_1& \cdots & {\widehat{T}}_{\delta }& 0_n& \cdots & 0_n& I_n& 0_n& \cdots & 0_n\\ 0_n& {\widehat{T}}_1& \cdots & {\widehat{T}}_{\delta }& 0_n& \cdots & 0_n& 0_n& 0_n& \cdots & 0_n\end{array}\right].$$

(19)

Remark 3.

Observe that increasing $n_b,$ with $n_b\ge \delta +1,$ the structures of W and ${\nabla }_Z(W)$ do not change due to the block band structure of the matrix W. Consequently, the length of the corresponding generators changes but their structure remains the same since only $T_0,T_1,\dots ,T_{\delta }$ and $I_n$ are different from zero in the first block row.

The computation by the GSA of the R factor of T and of its inverse $R^{-1}$ is made by only using the matrix G rather than the matrix T. Its implementation is a straightforward block matrix extension of the GSA described in Section 2.

Remark 4.

By construction, the initial generator matrix $G_0$ has the first $\delta +1$ block rows and the block row $n_b+1$ different from zero. Therefore, the multiplication of $G_0$ by the J-orthogonal matrix $H_1$ does not modify the structure of the generator matrix.

Let $G_0=G.$ At each iteration i (for $i=1,\dots ,n_b,$), we start from the generator matrix $G_{i-1}$ having the blocks (of length n) $i,i+1,\dots ,i+\delta$ and $n_b+1,\dots ,n_b+i$ different form zero. We then look for a J-orthogonal matrix $H_i$ such that the product $H_i{G}_{i-1}$ has in position $(1:n,(i-1)n+1:in)$ and $(n+1:2n,(i-1)n+1:in)$ a nonsingular upper triangular and zero matrix, respectively.

Then, $G_i$ is obtained from $\left[\begin{array}{c}{\tilde{G}}_i^{(+)}\\ {\tilde{G}}_i^{(-)}\end{array}\right]\equiv H_i{G}_{i-1}$ by multiplying the first n columns with $Z,$ i.e.,

$$G_i=\left[\begin{array}{c}{\tilde{G}}_i^{(+)}{Z}^T\\ {\tilde{G}}_i^{(-)}\end{array}\right].$$

The computation of the J-orthogonal matrix $H_i$ at the ith iteration of the GSA can be constructed as a product of n Householder matrices ${\widehat{H}}_{i,j}$ and n hyperbolic rotations ${\widehat{Y}}_{i,j},j=1,\dots ,n.$

The multiplication by the Householder matrices ${\widehat{H}}_{i,j}$ modifies the last n columns of the generator matrix, annihilating the last n entries but the $(n+1)$st in the row $(i-1)n+j,$ $j=1,\dots ,n,$ while the multiplication by the hyperbolic rotations ${\widehat{Y}}_{i,j}$ acts on the columns i and $n+1,$ annihilating the entry in position $((i-1)n+j,n+1).$

G_iven ${\upsilon }_1,{\upsilon }_2\in \mathbb{R},|{\upsilon }_1|>|{\upsilon }_2|,$ a hyperbolic matrix $Y\in {\mathbb{R}}^{2\times 2}$ can be computed

$$Y=\left[\begin{array}{cc}c& -s\\ -s& c\end{array}\right],\mathrm{with}c=\frac{{\upsilon }_1}{\sqrt{{\upsilon }_1^2-{\upsilon }_2^2}},s=\frac{{\upsilon }_2}{\sqrt{{\upsilon }_1^2-{\upsilon }_2^2}},$$

such that $[{\upsilon }_1,{\upsilon }_2]Y=[\sqrt{{\upsilon }_1^2-{\upsilon }_2^2},0].$

The modification of the sparsity pattern of the generator matrix after the first and ith iteration of the GSA are displayed in Figure 1 and  Figure 2, respectively.

The reliability of the GSA strongly depends on the way the hyperbolic rotation is computed. In [4,5,24], it is proven that the GSA is weakly stable if the hyperbolic rotations are implemented in an appropriate manner [3,11,12,24].

Let

$$H_{i,j}=\left[\begin{array}{cc}{I}_n& \\ & {\widehat{H}}_{i,j}\end{array}\right]Y_{i,j}=\left[\begin{array}{ccccc}{I}_{j-1}& & & & \\ & \hfill c_j& & \hfill -s_j& \\ & & I_{n-j}& & \\ & \hfill -s_j& & \hfill c_j& \\ & & & & I_{n-1}\end{array}\right].$$

Then,

$$H_i=H_{i,1}{Y}_{i,1}\cdots H_{i,n-1}{Y}_{i,n-1}{H}_{i,n}{Y}_{i,n}.$$

As previously mentioned, GSA relies only on the knowledge of the generators of W rather than on the matrix $\widehat{T}$ itself. Its computation involves the product ${\widehat{T}}^T\widehat{T},$ which can be accomplished with ${\delta }^2{n}^3$ flops. The ith iteration of the GSA involves the multiplication of n Householder matrices of size n times a matrix of size $((i+\delta +1)n\times n).$ Therefore, since the cost of the multiplication by the hyperbolic rotation is negligible with respect to that of the multiplication by the Householder matrices, the computational cost at iteration i is $4n^3(\delta +i)$. Hence, the computational cost of GSA is $2n^3{n}_b(2\delta +n_b^2/2)$.

5.2. GSA for Computing the Right Null-Space of Semidefinite Block Toeplitz Matrices

As already mentioned in Section 5.1, the number of desired blocks $n_b$ of the matrix T in Equation (13) can be computed as described in [8]. For the sake of simplicity, in the considered examples, we choose $n_b$ large enough to compute the null-space of $T.$

The structure and the computation via the GSA of the R factor of the $QR$ factorization of the singular block Toeplitz matrix T with rank $\rho <n\le m,$ is considered in [23].

A modification of the GSA for computing the null-space of Toeplitz matrices is described in [25]. In this paper, we extend the latter results to compute the null-space of T by modifying GSA.

Without loss of generality, let us assume that the first $\widehat{n}-1$ columns of T are linear independent and suppose that the $\widehat{n}$th column linearly depends on the previous ones. Therefore, the first $\widehat{n}-1$ principal minors of $\widehat{T}$ are positive while the $\widehat{n}$th one is zero. Let $\widehat{T}=Q\Lambda Q^T$ be the spectral decomposition of $\widehat{T},$ with $Q=[{\mathbf{q}}_1,\dots ,{\mathbf{q}}_{\widehat{n}}]$ orthogonal and $\Lambda =\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_{\widehat{n}-1},{\lambda }_{\widehat{n}}),$ with

$${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_{\widehat{n}-1}>{\lambda }_{\widehat{n}}=0,$$

and let ${\widehat{T}}_{\epsilon }=Q{\Lambda }_{\epsilon }{Q}^T,$ with ${\Lambda }_{\epsilon }=\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_{\widehat{n}-1},{\epsilon }^2),$ with $\epsilon \in {\mathbb{R}}_{+}^{\ast }.$ Hence,

$${\widehat{T}}_{\epsilon }^{-1}=\frac{1}{{\epsilon }^2}\left(\sum _{i=1}^{\widehat{n}-1}\frac{{\epsilon }^2}{{\lambda }_i}{\mathbf{q}}_i{\mathbf{q}}_i^T+{\mathbf{q}}_{\widehat{n}}{\mathbf{q}}_{\widehat{n}}^T\right).$$

Let $R_{\epsilon }$ be the Cholesky factor of ${\widehat{T}}_{\epsilon },$ with $R_{\epsilon }$ upper triangular, i.e., ${\widehat{T}}_{\epsilon }=R_{\epsilon }^T{R}_{\epsilon }.$ Then,

$${\epsilon }^2{\widehat{T}}_{\epsilon }^{-1}{e}_{\widehat{n}}^{(2\widehat{n})}=\left(\sum _{i=1}^{\widehat{n}-1}\frac{{\epsilon }^2{\mathbf{q}}_i^T{e}_{\widehat{n}}^{(2\widehat{n})}}{{\lambda }_i}{\mathbf{q}}_i+({\mathbf{q}}_{\widehat{n}}^T{e}_{\widehat{n}}^{(2\widehat{n})}){\mathbf{q}}_{\widehat{n}}\right).$$

On the other hand,

$${\epsilon }^2{\widehat{T}}_{\epsilon }^{-1}{e}_{\widehat{n}}^{(2\widehat{n})}={\epsilon }^2{R}_{\epsilon }^{-1}{R}_{\epsilon }^{-T}{e}_{\widehat{n}}^{(2\widehat{n})}={\epsilon }^2{r}_{\widehat{n},\widehat{n}}^{-1}{R}_{\epsilon }^{-1}{e}_{\widehat{n}}^{(2\widehat{n})},$$

where $r_{\widehat{n},\widehat{n}}={e_{\widehat{n}}^{(2\widehat{n})}}^T{R}_{\epsilon }{e}_{\widehat{n}}^{(2\widehat{n})}.$ Hence, as $\epsilon \to 0^{+},$ the last column of $R_{\epsilon }^{-1}$ becomes closer and closer to a multiple of ${\mathbf{q}}_{\widehat{n}},$ the eigenvector corresponding to the 0 eigenvalue of $\widehat{T}.$

Therefore, given

$$W_{\epsilon }=\left[\begin{array}{cc}{\widehat{T}}_{\epsilon }& I_{\widehat{n}}\\ I_{\widehat{n}}& 0_{\widehat{n}}\end{array}\right],$$

we have that

$${\nabla }_Z(W_{\epsilon })={G^{(+)}}^T{G}^{(+)}-{G^{(-)}}^T{G}^{(-)}+{\epsilon }^2{e}_{\widehat{n}}^{(2\widehat{n})}{e_{\widehat{n}}^{(2\widehat{n})}}^T.$$

Let

$$G_{0,\epsilon }=\left[\begin{array}{c}{G}^{(+)}\\ G^{(-)}\\ \epsilon {e_{\widehat{n}}^{(2\widehat{n})}}^T\end{array}\right].$$

Define

$$J_{\epsilon }=\mathrm{diag}(\underset{\widehat{n}}{\underbrace{1,1,\dots ,1}},\underset{\widehat{n}}{\underbrace{-1,-1,\dots ,-1}},1).$$

Hence,

$${\nabla }_Z(W_{\epsilon })=G_{0,\epsilon }^T{J}_{\epsilon }{G}_{0,\epsilon }.$$

We observe that column $\widehat{n}+1$ of the generator matrix is not involved in the GSA until the very last iteration, since only its $\widehat{n}$th entry is different from $0.$ At the very last iteration, the hyperbolic rotation

$$Y=\left[\begin{array}{cc}\hfill c& \hfill -s\\ \hfill -s& \hfill c\end{array}\right],$$

with

$$c=\frac{\sqrt{{G^{(+)}}^2(n,\widehat{n})+{\epsilon }^2}}{\sqrt{{G^{(+)}}^2(n,\widehat{n})+{\epsilon }^2-{G^{(-)}}^2(1,\widehat{n})}},s=\frac{G^{(-)}(1,\widehat{n})}{\sqrt{{G^{(+)}}^2(n,\widehat{n})+{\epsilon }^2-{G^{(-)}}^2(1,\widehat{n})}}$$

is applied to the $\widehat{n}$th and $(\widehat{n}+1)$st rows of $G,$ i.e., to the nth row of $G^{(+)}$ and the first one of $G^{(-)}$. Since $\widehat{T}$ is singular, it turns out that $|G^{(+)}(n,\widehat{n})|=|G^{(-)}(1,\widehat{n})|$ (see [23,25]). Thus,

$$\begin{array}{ccc}\hfill Y& =& \frac{1}{\sqrt{{G^{(+)}}^2(n,\widehat{n})+{\epsilon }^2-{G^{(-)}}^2(1,\widehat{n})}}\left[\begin{array}{cc}\sqrt{{G^{(+)}}^2(n,\widehat{n})+{\epsilon }^2}& -G^{(-)}(1,\widehat{n})\\ -G^{(-)}(1,\widehat{n})& \sqrt{{G^{(+)}}^2(n,\widehat{n})+{\epsilon }^2}\end{array}\right]\hfill \\ \hfill & =& \frac{|G^{(+)}(n,\widehat{n})|}{\epsilon }\left[\begin{array}{cc}\sqrt{1+{\left(\frac{\epsilon }{G^{(+)}(n,\widehat{n})}\right)}^2}& -\theta \\ -\theta & \sqrt{1+{\left(\frac{\epsilon }{G^{(+)}(n,\widehat{n})}\right)}^2}\end{array}\right],\hfill \end{array}$$

where

$$\theta =\frac{G^{(-)}(1,\widehat{n})}{|G^{(+)}(n,\widehat{n})|}=\mathrm{sign}(G^{(-)}(1,\widehat{n})).$$

We observe that, as $\epsilon \to 0^{+},$

$$\frac{|G^{(+)}(n,\widehat{n})|}{\epsilon }\to \infty ,\sqrt{1+{\left(\frac{\epsilon }{G^{(+)}(n,\widehat{n})}\right)}^2}\to 1.$$

Since a vector of the right null-space of T is determined except for the multiplication by a constant, neglecting the term $|G^{(+)}(n,\widehat{n})|/\epsilon ,$ such a vector can be computed at the last iteration as the first column of the product

$$\left[\begin{array}{cc}1& -\theta \\ -\theta & 1\end{array}\right]\left[\begin{array}{c}{G}^{(+)}(n,\widehat{n}+1:2\widehat{n})\\ G^{(-)}(1,\widehat{n}+1:2\widehat{n})\end{array}\right].$$

When detecting a vector of the null-space as a linear combination of row n of $G^{(+)}$ and row one of $G^{(-)},$ the new generator matrix G for the GSA is obtained removing the latter columns from G [23,25].

The implementation of the modified GSA for computing the null-space of band block-Toeplitz matrices in Equation (13) is rather technical and can be obtained from the authors upon request.

The stability properties of the GSA have been studied in [4,5,24]. The proposed algorithm inherits the stability properties of the GSA, which means that it is weakly stable.

6. Numerical Examples

All the numerical experiments were carried out in matlab with machine precision $\epsilon \approx 2.22\times {10}^{-16}$. Example 1 concerns the computation of the rank of a Sylvester matrix, while Examples 2 and 3 concern the computation of the null-space of polynomial matrices.

Example 1.

Let $x_i,i=1,\dots ,12,$ $y_i,i=1,\dots ,15,$ and $z_i,i=1,\dots ,3,$ be random numbers generated by thematlabfunctionrandn. Let $w(x)$ and $y(x)$ be the two polynomials of degree 15 and 18, constructed by the matlab functionpoly, whose roots are, respectively, $x_i$ and $z_j,i=1,\dots ,12,j=1,\dots ,3,$ and $y_i$ and $z_j,i=1,\dots ,15,j=1,\dots ,3.$

The greatest common divisor of w and y has degree 3 and, therefore, the Sylvester matrix $S\in {\mathbb{R}}^{33\times 33}$ constructed from $w(x)$ anf $y(x)$ has rank $30.$ The diagonal entries of the R factor computed by the $GSA$ implementation described in Section 4 are displayed in Figure 3. Observe that the rank of the matrix can be retrieved by the number of entries of R above a certain tolerance. Moreover, it can be noticed that the first $m=15$ diagonal entries monotonically decrease.

Example 2.

As second example, we consider the computation of the coprime factorization of a transfer function matrix, described in [9,26]. The results obtained by the proposed GSA-based algorithm were compared with those obtained computing the null-space of the considered matrix by the functionsvdofmatlab.

Let $H(s)=N_r(s)D_r^{-1}(s)$ be the transfer function with

$$D_r(s)=\left[\begin{array}{cccc}1-s& 0& 0& 0\\ 0& 1-s& 0& 0\\ 0& -s& 1-s& 0\\ 0& 0& 0& 1-s\end{array}\right],N_r(s)=\left[\begin{array}{cccc}{s}^2& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& s& 0\\ 0& 0& 0& s\end{array}\right].$$

Let

$$M(s)=\left[\begin{array}{cc}{N}_r^T(s)& -D_r^T(s)\end{array}\right]=\left[\begin{array}{ccccccccc}{s}^2& 0& 0& 0& 0& s-1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& s-1& s& 0\\ 0& 0& 0& s& 0& 0& 0& s-1& 0\\ 0& 0& 0& 0& s& 0& 0& 0& s-1\end{array}\right].$$

As reported in [9,26], a minimal polynomial basis for the right null-space of $M(s)$ is

$$N(s)=\left[\begin{array}{ccc}1-s& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& {(1-s)}^2& 0\\ 0& 0& 1-s\\ s^2& 0& 0\\ 0& s^2& 0\\ 0& s-s^2& 0\\ 0& 0& s\end{array}\right].$$

Let us consider $n_b=3.$ Then, $T\in {\mathbb{R}}^{20\times 21}$ is the block-Toeplitz matrix constructed from $M(s)$ as described in Section 5.1. Let $\mathrm{rank}(M)=17$ and $U\Sigma V^T$ be the rank and the singular value decomposition of T computed bymatlab, respectively, and let us define $V_1=V(:,1:17)$ and $V_2=V(:,18:21)$ the matrices of the right singular vectors corresponding to the nonzero and zero singular values of $T,$ respectively. The modified GSA applied to T yields four vectors ${\mathbf{v}}_1,Z_{n_b}{\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{v}}_3\in {\mathbb{R}}^{21}$ belonging to the right null-space of $M(s),$ with $Z_{n_b}=\mathrm{diag}(\mathrm{ones}(14,1),-7).$ Let $X=[{\mathbf{v}}_1{Z}_{n_b}{\mathbf{v}}_1{\mathbf{v}}_2{\mathbf{v}}_3].$ In

Table 1. Relative norm of $T{V}_2$, relative norm of $TX$, norm of $V_1^T{V}_2$ and norm of $V_1^{T}X,$ for Example 2.

$\frac{\parallel \mathit{T}{\mathit{V}}_2{\parallel }_2}{{\parallel \mathit{T}\parallel }_2}$	$\frac{{\parallel \mathit{T}\mathit{X}\parallel }_2}{{\parallel \mathit{T}\parallel }_2}$	$\parallel {\mathit{V}}_1^{\mathit{T}}{\mathit{V}}_2{\parallel }_2$	$\parallel {\mathit{V}}_1^{\mathit{T}}{\mathit{X}\parallel }_2$
$4.23\times {10}^{-16}$	$2.89\times {10}^{-16}$	$9.66\times {10}^{-16}$	$2.59\times {10}^{-15}$

, the relative norm of $T{V}_2$, the relative norm of $TX$, the norm of $V_1^T{V}_2$ and the norm of $V_1^{T}X,$ are reported in Columns 1–4, respectively.

Such values show that the results provided bysvdofmatlaband by the algorithm based on a modification of GSA are comparable in terms of accuracy.

Example 3.

This example can be found in [9,26]. Let $H(s)=D_l^{-1}(s)N_r(s)$ be the transfer function with

$$D_l(s)={(s+2)}^2(s+3)\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],N_l(s)=\left[\begin{array}{cc}3s+8& 2s^2+6s+2\\ s^2+6s+2& 3s^2+7s+8\end{array}\right].$$

Let $M(s)=[D_l(s),-N_L(s)].$ A right coprime pair for $M(s)$ is given by

$$N_r=\left[\begin{array}{cc}3& 2\\ s+2& 3\end{array}\right],D_r=\left[\begin{array}{cc}{s}^2+3s+4& 2\\ 2& s+4\end{array}\right],$$

Let us choose $n_b=4.$ Then, $T\in {\mathbb{R}}^{14\times 16}$ is the block-Toeplitz matrix constructed from $M(s)$ as described in Section 5.1. Let $\mathrm{rank}(M)=11$ and $U\Sigma V^T$ be the rank and the singular value decomposition of T computed bymatlab, respectively, and let define $V_1=V(:,1:11)$ and $V_2=V(:,12:16)$ the matrices of the right singular vectors corresponding to the nonzero and zero singular values of $T,$ respectively. The modified GSA applied to T yields the vectors ${\mathbf{v}}_1,Z_{n_b}{\mathbf{v}}_1,Z_{n_b}^2{\mathbf{v}}_1,{\mathbf{v}}_2,Z_{n_b}{\mathbf{v}}_2,{\mathbf{v}}_3\in {\mathbb{R}}^{16}$ of the right null-space, with $Z_{n_b}=\mathrm{diag}(\mathrm{ones}(12,1),-4).$ Let $X=[{\mathbf{v}}_1,Z_{n_b}{\mathbf{v}}_1,Z_{n_b}^2{\mathbf{v}}_1,{\mathbf{v}}_2,Z_{n_b}{\mathbf{v}}_2,{\mathbf{v}}_3].$ In

Table 2. Relative norm of $T{V}_2$, relative norm of $TX$, norm of $V_1^T{V}_2$ and norm of $V_1^{T}X,$ for Example 3.

$\frac{\parallel \mathit{T}{\mathit{V}}_2{\parallel }_2}{{\parallel \mathit{T}\parallel }_2}$	$\frac{{\parallel \mathit{T}\mathit{X}\parallel }_2}{{\parallel \mathit{T}\parallel }_2}$	$\parallel {\mathit{V}}_1^{\mathit{T}}{\mathit{V}}_2{\parallel }_2$	$\parallel {\mathit{V}}_1^{\mathit{T}}{\mathit{X}\parallel }_2$
$1.33\times {10}^{-16}$	$2.04\times {10}^{-16}$	$5.56\times {10}^{-16}$	$4.91\times {10}^{-15}$

, the relative norm of $T{V}_2$, the relative norm of $TX$, the norm of $V_1^T{V}_2$ and the norm of $V_1^{T}X,$ are reported in Columns 1–4, respectively.

As in Example 2, the results yielded by the considered algorithms are comparable in accuracy.

7. Conclusions

The Generalized Schur Algorithm is a powerful tool allowing to compute classical decompositions of matrices, such as the $QR$ and $LU$ factorizations. If the involved matrices have a particular structure, such as Toeplitz or Sylvester, the GSA computes the latter factorizations with a complexity of one order of magnitude less than that of classical algorithms based on Householder or elementary transformations.

After having emphasized the main features of the GSA, we have shown in this manuscript that the GSA helps to prove some theoretical properties of the R factor of the $QR$ factorization of some structured matrices. Moreover, a fast implementation of the GSA for computing the rank of Sylvester matrices and the null-space of polynomial matrices is proposed, which relies on a modification of the GSA for computing the R factor and its inverse of the $QR$ factorization of band block-Toeplitz matrices with full column rank. The numerical examples show that the proposed approach yields reliable results comparable to those ones provided by the function svd of matlab.