# Computing the Pseudoinverse of Specific Toeplitz Matrices Using Rank-One Updates

Authors: Predrag S. Stanimirović; Vasilios N. Katsikis; Igor Stojanović
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Application of the pure rank-one update algorithm as well as a combination of rank-one updates and the Sherman-Morrison formula in computing the Moore-Penrose inverse of the particular Toeplitz matrix is investigated in the present paper. Such Toeplitz matrices appear in the image restoration process and in many scientific areas that use the convolution. Four different approaches are developed, implemented, and tested on a number of numerical experiments.


---

## Body

## 1. Introduction

LetC
m
×
n and C
r
m
×
n denote the set of all complex m
×
n matrices and the set of all complex m
×
n matrices of rank
⁡
r, respectively. The identity matrix of an appropriate order is denoted by I. The conjugate transpose, the range, the rank, and the null space of A
∈
C
m
×
n are denoted by A
∗, R
(
A
), rank
⁡
A, and N
(
A
), respectively.Representation and computation of various generalized inverses are closely related to the following Penrose equations:(1)
A
X
A
=
A
,
X
A
X
=
X
,
A
X
∗
=
A
X
,
X
A
∗
=
X
A
. The set of all matrices obeying the conditions contained in a subset S
⊆
{
1,2
,
3,4
} is denoted by A
{
S
}. Any matrix from A
{
S
} is called S-inverse of A and is denoted by A
(
S
). By A
{
S
}
s we denote the set of all S-inverses of A with rank
⁡
s. For any matrix A there exists a single element in the set A
{
1,2
,
3,4
}, called the Moore-Penrose inverse of A and denoted by A
†.A rank-one modification of a matrixA
∈
C
m
×
n is the matrix M
=
A
+
b
c
∗, which is created from A and two vectors b
∈
C
m and c
∈
C
n. The Sherman-Morrison formula (S-M shortly) gives the basic relationship between the inverses M
-
1 and A
-
1 (for more details see, e.g., [1]):(2)
M
-
1
=
A
-
1
-
1
1
+
c
∗
A
-
1
b
A
-
1
b
c
∗
A
-
1
.The identity (2) provides a numerically efficient way to compute the inverse M
-
1 of the rank-one update M. The S-M formula is important in many different fields of numerical computation; see for example [2–7].On the other hand, Toeplitz matrices arise in a number of various theoretical investigations and applications. A number of iterative processes for finding generalized inverses of an arbitrary Toeplitz matrix by modifying Newton’s method have been developed so far. The main results were stated in [8–13]. Adaptations of the iterative processes to the Toeplitz structure are based on the usage of the displacement operator as well as the concept of displacement representation and ε-displacement rank of matrices.A variety of methods for computing the Moore-Penrose inverse of a rank-one modified matrix have been developed so far. Main results were derived in [14–18]. Relationships between various generalized inverses of an arbitrary matrix and corresponding generalized inverses of its rank-one modifications were investigated in [19]. The leading idea in [15, 16] was successive computation of the symmetric rank-one (SR1) updates A
i
-
1
+
a
i
a
i
∗
† of a given matrix A, where a
i
∗ denotes the ith row of A and A
i
-
1
=
∑
j
=
1
i
-
1
a
i
a
i
∗. The authors of the paper [15] introduced a computational procedure for the Moore-Penrose inverse of a symmetric rank-one perturbed matrix. Using this method, the authors of [16] proposed a finite method for computing the minimum-norm least-squares solution of the linear system A
x
=
b.The results derived in [20] reveal that both the SR1 updates techniques and the S-M recursive rule are useful tools in the computation of various matrix products involving the Moore-Penrose inverse of certain symmetric matrices. Particularly, the algorithms introduced in [20] are numerically efficient in computation of {
2,4
} and {
2,3
} inverses.In the present paper, we investigate the possibilities to apply the SR1 update procedure and the S-M formula in the computation of the Moore-Penrose inverse of specific Toeplitz matrices that appear in the image restoration process. Our main motivation arises from the convenience to apply the SR1 update and the S-M procedure in removing the blur which is always present in digital images. Firstly, both the SR1 update and the S-M formula are based on the usage of columns (or rows) of the input matrix. On the other hand, the matrices which appear in the mathematical model of blur in computer-generated images possess very specific structure which can be used to accelerate SR1 and S-M procedures. Namely, entries in Toeplitz matrices are constant along main diagonal parallels and, moreover, possess a significant proportion of zero elements.The paper is organized as follows. Some basic notations and necessary facts are restated in Section2. Also, some additional motivation is presented in the same section. Usage of the pure SR1 update algorithm, proposed in [15], in the computation of the Moore-Penrose inverse of a kind of Toeplitz matrices is considered in Section 3. A hybrid combination of the SR1 and the S-M recursive rules is defined in Section 4. An improvement of the SR1 procedure, which is derived on the basis of the specific structure of the underlying Toeplitz matrix, is presented in the same section. An application of introduced methods in image restoration is presented in Section 5.


## 2. Preliminaries and Motivation

Toeplitz matrices or diagonally constant matrices are matrices having constant diagonal entries. Toeplitz matrices which are applicable in the image restoration process containl nonzero main diagonal parallels above the main diagonal, where l defines the blurring process. In what follows, let us consider the Toeplitz matrix of such form:(3)
H
=
t
1
t
2
t
3
⋯
t
l
0
0
0
0
0
0
0
0
t
1
t
2
t
3
⋯
t
l
0
0
0
0
0
0
⋮
⋱
⋱
⋱
⋱
⋮
⋱
⋱
⋱
⋱
⋮
⋮
0
0
0
t
1
t
2
t
3
⋯
t
l
0
0
0
0
0
0
0
0
t
1
t
2
t
3
⋯
t
l
0
0
0
0
0
0
0
0
t
1
t
2
t
3
⋯
t
l
0
0
⋮
⋱
⋱
⋱
⋱
⋱
⋱
⋱
t
3
⋮
⋱
⋮
0
0
0
0
0
0
0
t
1
t
2
t
3
⋯
t
l
.The assumption t
1
≠
0 is active.To clarify notation, Toeplitz matrices of the general form (3) will be denoted shortly by (4)
T
p
l
t
1
,
t
1
,
…
,
t
l
.We investigate the use of the SR1 update method, as described in [15, Algorithm  2], during the numerical computation of the Moore-Penrose inverse of Toeplitz matrices satisfying the T
p
l
{
t
1
}
,
{
t
1
,
…
,
t
l
} pattern. Also, we examine different improvements of the original method. The improvements are based on appropriate adaptations of the SR1 method and the S-M formula to the characteristic structure of underlying matrices of type T
p
l
{
t
1
}
,
{
t
1
,
…
,
t
l
}. The method of SR1 updates is based on the expression which computes the Moore-Penrose inverse of the first k columns of the initial matrix A using the Moore-Penrose inverse of its first k
-
1 columns. In detail, the SR1 method from [15] starts from the well-known representation A
†
=
A
∗
A
†
A
∗ of the Moore-Penrose. If the ith row of A is denoted by a
i
∗, then(5)
A
∗
A
=
∑
i
=
1
m
a
i
a
i
∗
. Chen and Ji in [15] defined the matrix sequences A
k and X
k, as(6)
A
0
=
O
∈
C
n
×
n
,
A
k
=
∑
i
=
1
k
a
i
a
i
∗
,
k
=
1
,
…
,
m
,
X
k
=
A
k
†
A
∗
,
k
=
1
,
…
,
m
.Clearly, A
k
=
A
k
-
1
+
a
i
a
i
∗ is the rank-one modification of A
k
-
1 and(7)
X
m
=
A
m
†
A
∗
=
∑
i
=
1
m
a
i
a
i
∗
†
A
∗
=
A
†
.Recall that the Moore-Penrose inverse of a general rank-one modified matrixM
=
A
+
b
c
∗, where A is an arbitrary matrix and b,  c are arbitrary vectors, is obtained in [14, Theorem  3.1.3]. The general theorem from [14] suggests six different cases that one has to follow in order to establish a relation between M
† and A
†. In [15], the authors proved that the real number β
l
=
1
+
a
l
∗
A
l
-
1
†
a
l, corresponding to the term 1
+
c
∗
A
-
1
b in (2), satisfies β
l
≥
1. Later using this result in conjunction with the fact that A
l is a positive semidefinite matrix, the six cases of Theorem  3.1.3 from [14] can be reduced to the two-case problem. This reduction simplifies the SR1 updates formulas.Let us denote the firstk columns of H by(8)
H
k
=
h
1
⋯
h
k
. Also, the last n
-
k columns of H are denoted by {
k
+
1
}
H
{
n
}. Then the matrix H is given in the block form(9)
H
=
H
m
∣
H
n
m
+
1
∈
C
m
×
n
,
n
=
m
+
l
-
1
,
l
≥
1
,
H
m
∈
C
m
×
m
,
H
n
m
+
1
∈
C
m
×
n
-
m
, where the square block(10)
H
m
=
t
1
t
2
t
3
⋯
t
l
0
0
0
0
t
1
t
2
t
3
⋯
t
l
0
0
⋮
⋮
⋱
⋱
⋱
⋮
⋱
⋮
0
0
0
t
1
t
2
t
3
⋯
t
l
0
0
0
0
t
1
t
2
t
3
⋯
0
0
0
0
0
t
1
t
2
t
3
⋮
⋱
⋱
⋱
⋱
⋱
⋱
⋮
0
0
0
0
0
0
⋯
t
1 is the nonsingular band Toeplitz matrix and(11)
m
+
1
H
n
=
0
⋯
0
0
⋮
⋱
⋱
⋱
t
l
⋯
0
0
⋯
t
l
0
0
t
3
⋯
⋱
⋮
t
2
t
3
⋯
t
l collects the last n
-
m columns of H.An application of Greville’s partitioning method from [21] and the block partitioning method (BP method, shortly) from [22] in computing the Moore-Penrose inverse of the matrix H is presented in [23]. According to Algorithms  1 and 2 and Lemma  3 from [23], it is clear that the specific structure of matrix H enables computation H
† simply by computing the inverse of the nonsingular block H
{
m
}.In this paper, we investigate some alternative methods for computingH
† using the SR1 updates and the S-M formula. Our intention is to decrease computational complexity as much as possible using a specific structure of Toeplitz matrices of the general form T
p
l
{
t
1
}
,
{
t
1
,
…
,
t
l
}.


## 3. Computing the Pseudoinverse of a Toeplitz Matrix by Rank-One Updates

Our first attempt consists in applying the unmodified SR1 method from [15] in order to compute the pseudoinverse of the matrix H that belongs to the class T
p
l
{
t
1
}
,
{
t
1
,
…
,
t
l
}. Obtained results are compared with corresponding results derived by applying Algorithm  2 from [23] (the BP method). The results of this comparison are presented in Table 1, where n, l, and s are the parameters which define the Gaussian blur modeled by the Toeplitz matrix H. In the rest of the paper, it is assumed that the matrix H is of the order m
×
n, where n
=
m
+
l
-
1 and l
≥
1 represents the width of the blurring function.Table 1
Comparison of the rank-one updates method and the BP method.
Method
Sizesn, ℓ,  s
CPU time
r
1
r
2
r
3
r
4
BP
50, 15, 10
0.0025
2.0011e − 15
2.0385e − 14
3.773e − 15
1.6463e − 14
SR1
50, 15, 10
0.0490
2.2762e − 14
5.5984e − 13
1.2594e − 14
4.0623e − 13
BP
50, 20, 10
0.0028
6.9626e − 16
2.8774e − 14
3.4226e − 15
2.7706e − 14
SR1
50, 20, 10
0.0364
6.16e − 14
2.0705e − 12
1.3341e − 14
1.3718e − 12
BP
50, 15, 500
0.0029
8.0717e − 16
1.6354e − 14
3.7149e − 15
7.0839e − 15
SR1
50, 15, 500
0.0501
1.9666e − 14
4.8611e − 13
8.979e − 15
3.5898e − 13
BP
50, 20, 500
0.0016
6.8218e − 16
2.1407e − 14
3.2663e − 15
6.0107e − 15
SR1
50, 20, 500
0.0362
2.0844e − 14
5.3265e − 13
7.9594e − 15
4.1980e − 13
BP
250, 15, 500
0.0059
3.4724e − 15
2.3693e − 13
1.5467e − 14
7.4894e − 14
SR1
250, 15, 500
23.9762
7.81e − 12
3.6574e − 11
1.2707e − 13
2.2523e − 11
BP
400, 20, 500
0.0135
4.731e − 15
3.7769e − 13
2.6494e − 14
9.909e − 14
SR1
400, 20, 500
320.0578
3.8409e − 11
1.4246e − 10
2.2311e − 13
8.0923e − 11If an approximation ofH
† is denoted by X, then the residual norms are denoted by(12)
r
1
=
H
X
H
-
H
2
,
r
2
=
X
H
X
-
X
2
,
r
3
=
H
X
∗
-
H
X
2
,
r
4
=
X
H
∗
-
X
H
2
.For the sake of simplicity, let us denote Algorithm  2 from [15] by SR1 and Algorithm  2 from [23] by BP.From Table1, it is observable that the SR1 method is marginally efficient in terms of accuracy, but it is time consuming, especially for larger dimensions. There is a theoretical explanation for these results. The number of multiplications and divisions included in the SR1 algorithm is about O
(
m
2
n
+
m
n
) in the case m
<
n and O
(
n
2
m
+
m
n
) in the case n
<
m. Therefore, as it is stated in [15], the rank-one updates algorithm works efficiently in cases where it holds m
≪
n or m
≫
n. In our case, n
=
m
+
l
-
1, where l
≥
1 represents the width of the blurring process. Since, in the general case, l is a relatively small integer, the conditions m
<
n and m
≈
n are satisfied. This fact causes relative inefficiency of the SR1 method.Our second attempt is to apply Algorithm SR1 to the matrixH
m, from (10), in order to generate its inverse. Then the Moore-Penrose inverse of H comes from the block partitioning method (called BP), as it was described in [23]. This approach is shortly denoted by HSR1 + BP.It is observable from Table2 that the method HSR1 + BP is marginally efficient in terms of accuracy, but it is time consuming; see, for example, the case n
=
400.Table 2
Comparison of the combination of AlgSR1 and AlgBP methods and the block partitioning method.
Method
Sizesn, ℓ,  s
CPU time
r
1
r
2
r
3
r
4
BP
50, 15, 10
0.0037
2.0011e − 15
2.0385e − 14
3.773e − 15
1.6463e − 14
HSR1 + BP
50, 15, 10
0.0592
2.6951e − 14
2.9305e − 13
1.3171e − 13
4.0769e − 14
BP
50, 20, 10
3.4243e − 04
6.9626e − 16
2.8774e − 14
3.4226e − 15
2.7706e − 14
HSR1 + BP
50, 20, 10
0.0305
5.3084e − 14
4.1354e − 13
3.1072e − 13
9.9898e − 14
BP
50, 15, 500
0.0019
8.0717e − 16
1.6354e − 14
3.7149e − 15
7.0839e − 15
HSR1 + BP
50, 15, 500
0.0411
2.1453e − 14
2.2518e − 13
7.0401e − 14
3.2624e − 14
BP
50, 20, 500
3.2256e − 04
6.8218e − 16
2.1407e − 14
3.2663e − 15
6.0107e − 15
HSR1 + BP
50, 20, 500
0.0288
6.029e − 15
1.3385e − 13
2.6565e − 14
1.3927e − 14
BP
250, 15, 500
0.0048
3.4724e − 15
2.3693e − 13
1.5467e − 14
7.4894e − 14
HSR1 + BP
250, 15, 500
21.9978
1.8868e − 12
3.4807e − 10
1.154e − 11
4.0825e − 12
BP
400, 20, 500
0.0132
4.731e − 15
3.7769e − 13
2.6494e − 14
9.909e − 14
HSR1 + BP
400, 20, 500
306.1949
7.3251e − 12
2.7048e − 09
6.4634e − 11
2.0001e − 11From the previous analysis of data included in Tables1 and 2, we conclude that both algorithms SR1 and HSR1 + BP are not efficient approaches for computing H
† with respect to the BP method in terms of computational speed, especially for large matrices. But both approaches show marginal efficiency in terms of the accuracy.


## 4. A Combination of SR1 Updates and S-M Formula

Two additional algorithms are defined in the current section. The first one uses the SR1 updates to compute the matrixH
m
=
H
{
m
}
H
{
m
}
∗ in the first step and the S-M formula to compute H
m
† in the subsequent step. The second algorithm replaces the usage of the SR1 updates from the first step by a direct construction of the matrix H
m. The matrix H
m can be generated efficiently in view of the specific structure of the input matrix H.


### 4.1. SR1 Updates in Conjunction with the S-M Formula

As usual, theith column of H is denoted by h
i. Then, since H is of full row rank, it is advantageous to define a proper SR1 recursive computational method for generating the Moore-Penrose inverse of H on the basis of the relation(13)
H
†
=
H
∗
H
H
∗
-
1
=
H
∗
∑
i
=
1
n
h
i
h
i
∗
-
1
.For that purpose, it is necessary to define the matrix sequence(14)
H
0
=
O
∈
C
m
×
m
,
H
k
=
∑
i
=
1
k
h
i
h
i
∗
=
H
k
-
1
+
h
k
h
k
∗
,
k
=
1
,
…
,
n
.Lemma 1 reveals some of basic properties of the sequence H
k. Following the notation from [20], the notation h
k
∈
L
(
h
1
,
…
,
h
k
-
1
) means that a column h
k is linearly dependent on the previous columns h
1
,
…
,
h
k
-
1 and h
k
∉
L
(
h
1
,
…
,
h
k
-
1
) indicates that h
k is linearly independent of h
1
,
…
,
h
k
-
1.Lemma 1.
The following statements are valid:(i)
h
k
∉
L
(
h
1
,
…
,
h
k
-
1
), for each 1
<
k
≤
n.
(ii)
rank
⁡
H
k
=
rank
⁡
H
k
-
1
+
1, for each 1
<
k
≤
n.
(iii)
rank
⁡
H
k
<
m, for each 1
≤
k
<
m, and rank
⁡
H
k
=
rank
⁡
H
k
-
1
=
m, for each k
=
m
+
1
,
…
,
n.
(iv)
The following concluding relation is valid:(15)
H
∗
H
n
-
1
=
H
†
.Proof.
(i) This part of the proof is implied by the fact that the columns ofH are linearly independent.
(ii) According to the part (i), immediately followsrank
⁡
H
{
k
}
=
rank
⁡
H
{
k
-
1
}
+
1, for each 1
<
k
≤
n. Later, it is easy to verify(16)
H
k
=
∑
i
=
1
k
h
i
h
i
∗
=
h
1
⋯
h
k
h
1
∗
⋮
h
k
∗
=
H
k
H
k
∗
, which implies(17)
rank
⁡
H
k
=
rank
⁡
H
k
=
rank
⁡
H
k
-
1
+
1
=
rank
⁡
H
k
-
1
H
k
-
1
∗
+
1
=
rank
⁡
H
k
-
1
+
1
.
(iii) This part of the proof can be verified by using (i) and (ii), together with the fact thatH
k are all m
×
m matrices, for each k
=
1
,
…
,
n.
(iv) According to part (iii), matricesH
k, m
≤
k
≤
n, are invertible. Since(18)
H
n
=
H
n
H
n
∗
=
H
H
∗ is invertible, then it holds(19)
H
∗
H
n
-
1
=
H
∗
H
H
∗
-
1
=
H
†
, which completes the proof.Remark 2.
According to Lemma1 the matrices H
k are of the order m
×
m and it holds rank
⁡
H
k
<
m, for each 1
≤
k
<
m, while rank
⁡
H
k
=
rank
⁡
H
k
-
1
=
m, for each k
=
m
+
1
,
…
,
n. Since rank
⁡
H
n
=
rank
⁡
H
H
∗
=
m, the invertibility of H
n is guaranteed, so that we are able to apply the S-M formula which relates H
k
-
1 with H
k
-
1
-
1, for each k
=
m
+
1
,
…
,
n.
According to Lemma1, matrices H
k, m
≤
k
≤
n, are invertible. Therefore, it is possible to define the following matrix sequence:(20)
X
k
=
H
∗
H
k
-
1
,
k
=
m
,
…
,
n
.Identity (15) immediately implies X
n
=
H
†
.Lemma 3.
It holds1
+
h
k
∗
H
k
-
1
-
1
h
k
≥
1, for each k
=
m
+
1
,
…
,
n.Proof.
The proof of lemma follows from the fact thatH
k
-
1 is regular, for each k
=
m
+
1
,
…
,
n.Now, we turn our attention to (20). Since the matrices H
k
-
1, k
=
m
+
1
,
…
,
n are regular, it is possible to find a relation between H
k
-
1
=
H
k
-
1
+
h
k
h
k
∗
-
1 and H
k
-
1
-
1 by applying the Sherman-Morrison formula (2). This relation is given in Proposition 4.Proposition 4.
Let the matrixH be defined as in (3). Assume that the matrix sequence H
k is defined in (14). The inverse of H
k is equal to(21)
H
k
-
1
=
H
k
-
1
+
h
k
h
k
∗
-
1
=
H
k
-
1
-
1
-
1
+
h
k
∗
H
k
-
1
-
1
h
k
-
1
H
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
,for each k
=
m
+
1
,
…
,
n.In order to present the finite recursive algorithm for computing the Moore-Penrose inverse ofH with the help of the actual SR1 update formulas (20) and (21), we first define H
0
=
O
∈
C
m
×
m. For each k
=
1
,
…
,
m, it is necessary to compute H
k
=
H
k
-
1
+
h
k
h
k
∗. This requires calculated values of the vectors y
m
,
t
=
H
m
-
1
h
t
∈
C
m, t
=
1
,
…
,
n and the matrix X
m
=
H
∗
H
m
-
1.Theorem 5.
Let the matrixH be defined as in (3), the matrix sequence H
k defined in (14), and the matrix sequence X
k defined in (20). Then the following recursive relations are true:(22)
y
k
,
t
=
y
k
-
1
,
t
-
y
k
-
1
,
k
h
k
∗
y
k
-
1
,
t
1
+
h
k
∗
y
k
-
1
,
k
=
I
m
-
y
k
-
1
,
k
h
k
∗
1
+
h
k
∗
y
k
-
1
,
k
·
y
k
-
1
,
t
,
(23)
X
k
=
X
k
-
1
-
H
∗
y
k
-
1
,
k
y
k
-
1
,
k
∗
1
+
h
k
∗
y
k
-
1
,
k
.Proof.
Note that fromy
l
,
t
=
H
l
-
1
h
t it follows that (
y
l
,
t
)
∗
=
H
l
-
1
h
t
∗
=
h
t
∗
H
l
-
1. Based on the use of (21), the recursive sequences y
k
,
t are defined for k
=
m
+
1
,
…
,
n as in the following:(24)
y
k
,
t
=
H
k
-
1
h
t
=
H
k
-
1
+
h
k
h
k
∗
-
1
h
t
=
H
k
-
1
-
1
-
1
+
h
k
∗
H
k
-
1
-
1
h
k
-
1
H
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
h
t
=
H
k
-
1
-
1
h
t
-
H
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
h
t
1
+
h
k
∗
H
k
-
1
-
1
h
k
,
t
=
1
,
…
,
n
.Then the statement (22) can be easily verified using (24).
Also, matricesX
k, k
=
m
+
1
,
…
,
n are defined as (25)
X
k
=
H
∗
H
k
-
1
=
H
∗
H
k
-
1
-
1
-
1
+
h
k
∗
H
k
-
1
-
1
h
k
-
1
H
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
=
H
∗
H
k
-
1
-
1
-
H
∗
1
+
h
k
∗
H
k
-
1
-
1
h
k
-
1
G
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
=
X
k
-
1
-
H
∗
H
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
1
+
h
k
∗
H
k
-
1
-
1
h
k
.So, (23) holds.According to Theorem5, we present Algorithm 1 for computing the Moore-Penrose inverse of the Toeplitz matrix H. The algorithm is based on the SR1 modifications to compute H
m and subsequently the S-M formula to compute H
†. For this purpose, we use the abbreviation HSR1 + S-M to denote Algorithm 1.Algorithm 1:Computing the Moore-Penrose inverse of the Toeplitz matrix H.
Require: Start with H
∈
C
m
×
n and n
>
m. Let h
i be the ith column of H, for i
=
1
,
…
,
n, H
0
=
O
∈
C
m
×
m.
(1) (SR1 Step) Fork
=
1
,
…
,
m, compute H
k
=
H
k
-
1
+
h
k
h
k
∗.
(2) Sety
m
,
t
=
H
m
-
1
h
t
∈
C
m for all t
=
1
,
…
,
n and X
m
=
H
∗
H
m
-
1.
(3) (S-M Step) Computey
l
,
t and X
l by using (22) and (23), respectively for all l
=
m
+
1
,
…
,
n.
(4) OutputH
†
=
X
n.AMatlab implementation of Algorithm 1 is placed in Appendix.


### 4.2. An Improvement of the Hybrid Combination in Algorithm1

By taking into account the fact that the Toeplitz matrixH has a specific structure, the matrix H
m can be generated in an efficient way. In this way, an improvement of step  (1) of Algorithm 1 is introduced. The effectiveness of this method will be justified in the numerical experiments presented in Section 4.3.According to the proper structure of the matrixH and step  (1) of Algorithm 1, the general structure of the matrix H
m is defined as(26)Let w
=
(
w
1
,
…
,
w
n
) be one-dimensional vector; then the symmetric Toeplitz matrix generated by w is denoted by T
o
e
p
l
i
t
z
(
w
):(27)
T
o
e
p
l
i
t
z
w
=
w
1
w
2
w
3
⋯
w
n
w
2
w
1
w
2
⋯
w
n
-
1
⋮
⋱
⋱
⋱
⋮
w
n
-
1
⋯
⋯
w
1
w
2
w
n
⋯
⋯
w
2
w
1
.Also, for two appropriate vectors u
,
v, we denote by F
=
H
a
n
k
e
l
(
u
,
v
) the (symmetric) Hankel matrix whose first column is u and whose last row is v.Theorem 6.
The matrixH
m can be expressed as the difference of a specific banded symmetric Toeplitz matrix and a partitioned matrix of the form(28)
O
m
-
l
+
1
O
m
-
l
+
1
×
l
-
1
O
l
-
1
×
m
-
l
+
1
F
l
-
1
,where the block F
l
-
1 is defined by F
l
-
1
=
F
2, for an appropriate Hankel matrix F.Proof.
Letu and v be vectors defined by(29)
u
=
t
l
t
l
-
1
⋯
t
2
,
(30)
v
=
0
⋯
0
︸
l
-
2
t
l
,respectively. Then we define the following matrices:(31)
E
=
T
o
e
p
l
i
t
z
∑
i
=
1
l
t
i
2
∑
i
=
1
l
-
1
t
i
t
i
+
1
∑
i
=
1
l
-
2
t
i
t
i
+
2
⋯
t
1
t
l
0
⋯
0
︸
m
-
l
,
F
=
H
a
n
k
e
l
v
,
u
,
L
=
O
m
-
l
+
1
O
m
-
l
+
1
×
l
-
1
O
l
-
1
×
m
-
l
+
1
F
l
-
1
,
F
l
-
1
=
F
2
. A simple verification shows that H
m
=
E
-
L.Next we present an example, in low dimensions, in order to illustrate Theorem6.Example 7.
Consider the Toeplitz matrixH
∈
C
m
×
n, defined in (3) for m
=
5,  l
=
3, and n
=
m
+
l
-
1
=
7: (32)
H
=
H
5
∣
H
7
6
=
t
1
t
2
t
3
0
0
0
0
0
t
1
t
2
t
3
0
0
0
0
0
t
1
t
2
t
3
0
0
0
0
0
t
1
t
2
t
3
0
0
0
0
0
t
1
t
2
t
3
.
According to (14),(33)
H
m
=
H
5
=
H
5
H
5
∗
=
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
0
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
t
1
t
3
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
t
1
t
3
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
t
1
t
2
0
0
t
1
t
3
t
1
t
2
t
1
2
.We exploit the matrix(34)
E
=
T
o
e
p
l
i
t
z
∑
i
=
1
3
t
i
2
∑
i
=
1
2
t
i
t
i
+
1
t
1
t
3
0
0
=
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
0
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
t
1
t
3
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
t
1
t
3
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
0
0
t
1
t
3
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
.Letu
=
t
3
t
2 and v
=
0
t
3; then(35)
F
=
H
a
n
k
e
l
v
,
u
=
0
t
3
t
3
t
2
,
F
2
=
t
3
2
t
2
t
3
t
2
t
3
t
2
2
+
t
3
2
. Finally, it is necessary to construct the matrix(36)
L
=
O
3
O
3
×
2
O
2
×
3
F
2
=
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
t
3
2
t
2
t
3
0
0
0
t
2
t
3
t
2
2
+
t
3
2
.It is easy to see that H
5
=
E
-
L
.According to Theorem6 we present the following improved version of Algorithm 1, called Algorithm 2. AMatlab implementation of Algorithm 2 is placed in Appendix. It is appropriate to use the term IHSR1 to denote the improvement of the HSR1 step of Algorithm 1 and IHSR1 + S-M to denote Algorithm 2.Algorithm 2:Computing the Moore-Penrose inverse of the Toeplitz matrix H.
Require: Start with H
∈
C
m
×
n, n
>
m. Let h
i be the ith column of H, for i
=
1
,
…
,
n.
(1) (Improved HSR1 Step) DefineE and L as in Theorem 6. Then, set H
m
=
E
-
L.
(2) Sety
m
,
t
=
H
m
-
1
h
t
∈
C
m for all t
=
1
,
…
,
n and X
m
=
H
∗
H
m
-
1.
(3) (S-M Step) Computey
l
,
t and X
l by using (22) and (23), respectively for all l
=
m
+
1
,
…
,
n.
(4) OutputH
†
=
X
n.Note that proposed Algorithms1 and 2 are not using directly, at any step, the Sherman-Morrison formula to compute the inverse of some matrix. In Section 4.4, we provide more information regarding this situation.


### 4.3. Numerical Experiments

In this subsection we analyze numerical data arising during the computation of the Moore-Penrose inverse of the Toeplitz matrixH by applying aMatlab implementation of Algorithms 1 and 2. In order to test the time efficiency as well as the accuracy of considered methods, we enrich our collection of matrices used in Tables 1 and 2 with larger matrices. In order to complete our numerical study, we also compare the results derived by applying Algorithm  2 from [23] (Algorithm BP). A comparison of Algorithms 1 and 2 and BP is presented in Table 3.Table 3
Comparison of Algorithm1 (HSR1 + S-M), Algorithm 2 (IHSR1 + S-M), and BP.
Method
Sizesn, ℓ, s
CPU time
r
1
r
2
r
3
r
4
HSR1 + S-M
50, 15, 10
0.0054
1.6449e − 14
2.5845e − 12
3.5344e − 13
1.9362e − 14
IHSR1 + S-M
50, 15, 10
0.0035
2.0486e − 14
2.7318e − 12
5.4426e − 13
2.3222e − 14
BP
50, 15, 10
0.0018
2.0011e − 15
2.0385e − 14
3.773e − 15
1.6463e − 14
HSR1 + S-M
50, 20, 10
0.0040
1.8748e − 14
8.8674e − 13
2.2185e − 13
1.5379e − 14
IHSR1 + S-M
50, 20, 10
0.0032
9.4733e − 15
8.9539e − 13
1.9076e − 13
1.4775e − 14
BP
50, 20, 10
0.0027
9.2684e − 16
1.6646e − 14
3.2679e − 15
5.8648e − 15
HSR1 + S-M
50, 20, 500
0.0053
1.2862e − 14
9.6327e − 13
1.2288e − 13
1.4447e − 14
IHSR1 + S-M
50, 20, 500
0.0044
8.6971e − 15
1.5068e − 12
2.2331e − 13
1.5216e − 14
BP
50, 20, 500
0.0022
6.8218e − 16
2.1407e − 14
3.2663e − 15
6.0107e − 15
HSR1 + S-M
250, 15, 500
0.0482
2.5535e − 13
3.7054e − 11
4.7471e − 12
2.5608e − 13
IHSR1 + S-M
250, 15, 500
0.0195
2.209e − 13
3.7065e − 11
5.5401e − 12
2.4239e − 13
BP
250, 15, 500
0.0032
3.4724e − 15
2.3693e − 13
1.5467e − 14
7.4894e − 14
HSR1 + S-M
400, 20, 500
0.3234
4.5138e − 13
6.0769e − 10
1.4669e − 11
5.5243e − 13
IHSR1 + S-M
400, 20, 500
0.1466
5.881e − 13
6.0878e − 10
1.5866e − 11
5.3143e − 13
BP
400, 20, 500
0.0116
6.5483e − 15
6.052e − 13
3.0382e − 14
1.3117e − 13
HSR1 + S-M
1200, 20, 500
11.0062
3.0653e − 12
1.7107e − 08
1.4049e − 10
3.2612e − 12
IHSR1 + S-M
1200, 20, 500
2.3018
3.5064e − 12
1.711e − 08
1.3136e − 10
3.1795e − 12
BP
1200, 20, 500
0.1049
1.2805e − 14
3.179e − 12
8.9833e − 14
2.3321e − 13
HSR1 + S-M
1500, 20, 500
21.7539
4.3812e − 12
3.2457e − 08
2.2848e − 10
4.2712e − 12
IHSR1 + S-M
1500, 20, 500
3.7358
4.2345e − 12
3.2456e − 08
2.0156e − 10
4.2919e − 12
BP
1500, 20, 500
0.1886
1.3786e − 14
4.0422e − 12
1.1081e − 13
3.0133e − 13
HSR1 + S-M
2000, 50, 500
74.5230
5.5113e − 12
4.6497e − 08
5.1185e − 10
6.7183e − 12
IHSR1 + S-M
2000, 50, 500
30.8622
5.3682e − 12
4.6504e − 08
5.4329e − 10
6.6102e − 12
BP
2000, 50, 500
0.4038
1.8118e − 14
6.5462e − 12
1.408e − 13
3.6715e − 13From Table3, it is observable that the level of improvement in CPU time that can be achieved by using Algorithm 2 instead of Algorithm 1 is significant. Also, it is clear that the accuracy of both algorithms is almost the same. Clearly, the BP method requires minimal CPU time.


### 4.4. Rounding Error Analysis

The starting motivation of this subsection is the following basic problem: approximation error grows with repeated use of the S-M formula. As a consequence, computations which involve the S-M rule become unstable. For this purpose, it is interesting to investigate the error curves corresponding to the four Penrose equations during the iterations included in step  (3) of Algorithm2. Nevertheless, it is important to note that the proposed Algorithm 2 does not use in each step the S-M rule. Of course, Theorem 5 is in the heart of Algorithm 2 but also it is clear that the poor stability of the S-M formula has no serious influence on the stability of proposed algorithms.So, in this subsection, we shall present error estimations regarding Algorithm2. That is, we record and investigate the residual matrix norms r
1
,
r
2
,
r
3, and r
4, corresponding to the four Penrose equations during the execution of step  (3) of Algorithm 2. The results are presented in Figures 1, 2, and 3 for the cases: (a) n
=
1200, l
=
20, and s
=
500, (b) n
=
1500, l
=
20, and s
=
500, and (c) n
=
2000, l
=
20
,
a
n
d
s
=
500, respectively. The choice of these matrices from Table 3 is random.Figure 1
Errors recorded forn
=
1200, l
=
20, and m
=
500 while the loop of step (3) is iterating.Figure 2
Errors recorded forn
=
1500, l
=
20, and m
=
500 while the loop of step (3) is iterating.Figure 3
Errors recorded forn
=
2000, l
=
50, and m
=
500 while the loop of step (3) is iterating.The purple line in each of Figures1–3 corresponds to the values of the matrix norm (
X
l
H
)
∗
-
X
l
H
2, l
=
m
+
1
,
…
,
n. Since X
l
=
H
∗
H
l
-
1 and the matrix H
l is symmetric, the matrix X
l
H is symmetric. For this reason, the graph corresponding to these values is almost a constant line.The stable convergence is observable from each of these figures. Also, it is possible to notice very fast convergence in the terminal phase of the convergence as well as a relatively slower convergence in the middle of the recurrent process. This fact is understandable if one takes into account that it is necessary to take into consideration all the columns of the input matrix in order to get the complete information on the pseudoinverse.


## 4.1. SR1 Updates in Conjunction with the S-M Formula

As usual, theith column of H is denoted by h
i. Then, since H is of full row rank, it is advantageous to define a proper SR1 recursive computational method for generating the Moore-Penrose inverse of H on the basis of the relation(13)
H
†
=
H
∗
H
H
∗
-
1
=
H
∗
∑
i
=
1
n
h
i
h
i
∗
-
1
.For that purpose, it is necessary to define the matrix sequence(14)
H
0
=
O
∈
C
m
×
m
,
H
k
=
∑
i
=
1
k
h
i
h
i
∗
=
H
k
-
1
+
h
k
h
k
∗
,
k
=
1
,
…
,
n
.Lemma 1 reveals some of basic properties of the sequence H
k. Following the notation from [20], the notation h
k
∈
L
(
h
1
,
…
,
h
k
-
1
) means that a column h
k is linearly dependent on the previous columns h
1
,
…
,
h
k
-
1 and h
k
∉
L
(
h
1
,
…
,
h
k
-
1
) indicates that h
k is linearly independent of h
1
,
…
,
h
k
-
1.Lemma 1.
The following statements are valid:(i)
h
k
∉
L
(
h
1
,
…
,
h
k
-
1
), for each 1
<
k
≤
n.
(ii)
rank
⁡
H
k
=
rank
⁡
H
k
-
1
+
1, for each 1
<
k
≤
n.
(iii)
rank
⁡
H
k
<
m, for each 1
≤
k
<
m, and rank
⁡
H
k
=
rank
⁡
H
k
-
1
=
m, for each k
=
m
+
1
,
…
,
n.
(iv)
The following concluding relation is valid:(15)
H
∗
H
n
-
1
=
H
†
.Proof.
(i) This part of the proof is implied by the fact that the columns ofH are linearly independent.
(ii) According to the part (i), immediately followsrank
⁡
H
{
k
}
=
rank
⁡
H
{
k
-
1
}
+
1, for each 1
<
k
≤
n. Later, it is easy to verify(16)
H
k
=
∑
i
=
1
k
h
i
h
i
∗
=
h
1
⋯
h
k
h
1
∗
⋮
h
k
∗
=
H
k
H
k
∗
, which implies(17)
rank
⁡
H
k
=
rank
⁡
H
k
=
rank
⁡
H
k
-
1
+
1
=
rank
⁡
H
k
-
1
H
k
-
1
∗
+
1
=
rank
⁡
H
k
-
1
+
1
.
(iii) This part of the proof can be verified by using (i) and (ii), together with the fact thatH
k are all m
×
m matrices, for each k
=
1
,
…
,
n.
(iv) According to part (iii), matricesH
k, m
≤
k
≤
n, are invertible. Since(18)
H
n
=
H
n
H
n
∗
=
H
H
∗ is invertible, then it holds(19)
H
∗
H
n
-
1
=
H
∗
H
H
∗
-
1
=
H
†
, which completes the proof.Remark 2.
According to Lemma1 the matrices H
k are of the order m
×
m and it holds rank
⁡
H
k
<
m, for each 1
≤
k
<
m, while rank
⁡
H
k
=
rank
⁡
H
k
-
1
=
m, for each k
=
m
+
1
,
…
,
n. Since rank
⁡
H
n
=
rank
⁡
H
H
∗
=
m, the invertibility of H
n is guaranteed, so that we are able to apply the S-M formula which relates H
k
-
1 with H
k
-
1
-
1, for each k
=
m
+
1
,
…
,
n.
According to Lemma1, matrices H
k, m
≤
k
≤
n, are invertible. Therefore, it is possible to define the following matrix sequence:(20)
X
k
=
H
∗
H
k
-
1
,
k
=
m
,
…
,
n
.Identity (15) immediately implies X
n
=
H
†
.Lemma 3.
It holds1
+
h
k
∗
H
k
-
1
-
1
h
k
≥
1, for each k
=
m
+
1
,
…
,
n.Proof.
The proof of lemma follows from the fact thatH
k
-
1 is regular, for each k
=
m
+
1
,
…
,
n.Now, we turn our attention to (20). Since the matrices H
k
-
1, k
=
m
+
1
,
…
,
n are regular, it is possible to find a relation between H
k
-
1
=
H
k
-
1
+
h
k
h
k
∗
-
1 and H
k
-
1
-
1 by applying the Sherman-Morrison formula (2). This relation is given in Proposition 4.Proposition 4.
Let the matrixH be defined as in (3). Assume that the matrix sequence H
k is defined in (14). The inverse of H
k is equal to(21)
H
k
-
1
=
H
k
-
1
+
h
k
h
k
∗
-
1
=
H
k
-
1
-
1
-
1
+
h
k
∗
H
k
-
1
-
1
h
k
-
1
H
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
,for each k
=
m
+
1
,
…
,
n.In order to present the finite recursive algorithm for computing the Moore-Penrose inverse ofH with the help of the actual SR1 update formulas (20) and (21), we first define H
0
=
O
∈
C
m
×
m. For each k
=
1
,
…
,
m, it is necessary to compute H
k
=
H
k
-
1
+
h
k
h
k
∗. This requires calculated values of the vectors y
m
,
t
=
H
m
-
1
h
t
∈
C
m, t
=
1
,
…
,
n and the matrix X
m
=
H
∗
H
m
-
1.Theorem 5.
Let the matrixH be defined as in (3), the matrix sequence H
k defined in (14), and the matrix sequence X
k defined in (20). Then the following recursive relations are true:(22)
y
k
,
t
=
y
k
-
1
,
t
-
y
k
-
1
,
k
h
k
∗
y
k
-
1
,
t
1
+
h
k
∗
y
k
-
1
,
k
=
I
m
-
y
k
-
1
,
k
h
k
∗
1
+
h
k
∗
y
k
-
1
,
k
·
y
k
-
1
,
t
,
(23)
X
k
=
X
k
-
1
-
H
∗
y
k
-
1
,
k
y
k
-
1
,
k
∗
1
+
h
k
∗
y
k
-
1
,
k
.Proof.
Note that fromy
l
,
t
=
H
l
-
1
h
t it follows that (
y
l
,
t
)
∗
=
H
l
-
1
h
t
∗
=
h
t
∗
H
l
-
1. Based on the use of (21), the recursive sequences y
k
,
t are defined for k
=
m
+
1
,
…
,
n as in the following:(24)
y
k
,
t
=
H
k
-
1
h
t
=
H
k
-
1
+
h
k
h
k
∗
-
1
h
t
=
H
k
-
1
-
1
-
1
+
h
k
∗
H
k
-
1
-
1
h
k
-
1
H
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
h
t
=
H
k
-
1
-
1
h
t
-
H
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
h
t
1
+
h
k
∗
H
k
-
1
-
1
h
k
,
t
=
1
,
…
,
n
.Then the statement (22) can be easily verified using (24).
Also, matricesX
k, k
=
m
+
1
,
…
,
n are defined as (25)
X
k
=
H
∗
H
k
-
1
=
H
∗
H
k
-
1
-
1
-
1
+
h
k
∗
H
k
-
1
-
1
h
k
-
1
H
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
=
H
∗
H
k
-
1
-
1
-
H
∗
1
+
h
k
∗
H
k
-
1
-
1
h
k
-
1
G
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
=
X
k
-
1
-
H
∗
H
k
-
1
-
1
h
k
h
k
∗
H
k
-
1
-
1
1
+
h
k
∗
H
k
-
1
-
1
h
k
.So, (23) holds.According to Theorem5, we present Algorithm 1 for computing the Moore-Penrose inverse of the Toeplitz matrix H. The algorithm is based on the SR1 modifications to compute H
m and subsequently the S-M formula to compute H
†. For this purpose, we use the abbreviation HSR1 + S-M to denote Algorithm 1.Algorithm 1:Computing the Moore-Penrose inverse of the Toeplitz matrix H.
Require: Start with H
∈
C
m
×
n and n
>
m. Let h
i be the ith column of H, for i
=
1
,
…
,
n, H
0
=
O
∈
C
m
×
m.
(1) (SR1 Step) Fork
=
1
,
…
,
m, compute H
k
=
H
k
-
1
+
h
k
h
k
∗.
(2) Sety
m
,
t
=
H
m
-
1
h
t
∈
C
m for all t
=
1
,
…
,
n and X
m
=
H
∗
H
m
-
1.
(3) (S-M Step) Computey
l
,
t and X
l by using (22) and (23), respectively for all l
=
m
+
1
,
…
,
n.
(4) OutputH
†
=
X
n.AMatlab implementation of Algorithm 1 is placed in Appendix.


## 4.2. An Improvement of the Hybrid Combination in Algorithm1

By taking into account the fact that the Toeplitz matrixH has a specific structure, the matrix H
m can be generated in an efficient way. In this way, an improvement of step  (1) of Algorithm 1 is introduced. The effectiveness of this method will be justified in the numerical experiments presented in Section 4.3.According to the proper structure of the matrixH and step  (1) of Algorithm 1, the general structure of the matrix H
m is defined as(26)Let w
=
(
w
1
,
…
,
w
n
) be one-dimensional vector; then the symmetric Toeplitz matrix generated by w is denoted by T
o
e
p
l
i
t
z
(
w
):(27)
T
o
e
p
l
i
t
z
w
=
w
1
w
2
w
3
⋯
w
n
w
2
w
1
w
2
⋯
w
n
-
1
⋮
⋱
⋱
⋱
⋮
w
n
-
1
⋯
⋯
w
1
w
2
w
n
⋯
⋯
w
2
w
1
.Also, for two appropriate vectors u
,
v, we denote by F
=
H
a
n
k
e
l
(
u
,
v
) the (symmetric) Hankel matrix whose first column is u and whose last row is v.Theorem 6.
The matrixH
m can be expressed as the difference of a specific banded symmetric Toeplitz matrix and a partitioned matrix of the form(28)
O
m
-
l
+
1
O
m
-
l
+
1
×
l
-
1
O
l
-
1
×
m
-
l
+
1
F
l
-
1
,where the block F
l
-
1 is defined by F
l
-
1
=
F
2, for an appropriate Hankel matrix F.Proof.
Letu and v be vectors defined by(29)
u
=
t
l
t
l
-
1
⋯
t
2
,
(30)
v
=
0
⋯
0
︸
l
-
2
t
l
,respectively. Then we define the following matrices:(31)
E
=
T
o
e
p
l
i
t
z
∑
i
=
1
l
t
i
2
∑
i
=
1
l
-
1
t
i
t
i
+
1
∑
i
=
1
l
-
2
t
i
t
i
+
2
⋯
t
1
t
l
0
⋯
0
︸
m
-
l
,
F
=
H
a
n
k
e
l
v
,
u
,
L
=
O
m
-
l
+
1
O
m
-
l
+
1
×
l
-
1
O
l
-
1
×
m
-
l
+
1
F
l
-
1
,
F
l
-
1
=
F
2
. A simple verification shows that H
m
=
E
-
L.Next we present an example, in low dimensions, in order to illustrate Theorem6.Example 7.
Consider the Toeplitz matrixH
∈
C
m
×
n, defined in (3) for m
=
5,  l
=
3, and n
=
m
+
l
-
1
=
7: (32)
H
=
H
5
∣
H
7
6
=
t
1
t
2
t
3
0
0
0
0
0
t
1
t
2
t
3
0
0
0
0
0
t
1
t
2
t
3
0
0
0
0
0
t
1
t
2
t
3
0
0
0
0
0
t
1
t
2
t
3
.
According to (14),(33)
H
m
=
H
5
=
H
5
H
5
∗
=
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
0
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
t
1
t
3
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
t
1
t
3
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
t
1
t
2
0
0
t
1
t
3
t
1
t
2
t
1
2
.We exploit the matrix(34)
E
=
T
o
e
p
l
i
t
z
∑
i
=
1
3
t
i
2
∑
i
=
1
2
t
i
t
i
+
1
t
1
t
3
0
0
=
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
0
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
t
1
t
3
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
t
1
t
3
0
t
1
t
3
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
t
1
t
2
+
t
2
t
3
0
0
t
1
t
3
t
1
t
2
+
t
2
t
3
t
1
2
+
t
2
2
+
t
3
2
.Letu
=
t
3
t
2 and v
=
0
t
3; then(35)
F
=
H
a
n
k
e
l
v
,
u
=
0
t
3
t
3
t
2
,
F
2
=
t
3
2
t
2
t
3
t
2
t
3
t
2
2
+
t
3
2
. Finally, it is necessary to construct the matrix(36)
L
=
O
3
O
3
×
2
O
2
×
3
F
2
=
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
t
3
2
t
2
t
3
0
0
0
t
2
t
3
t
2
2
+
t
3
2
.It is easy to see that H
5
=
E
-
L
.According to Theorem6 we present the following improved version of Algorithm 1, called Algorithm 2. AMatlab implementation of Algorithm 2 is placed in Appendix. It is appropriate to use the term IHSR1 to denote the improvement of the HSR1 step of Algorithm 1 and IHSR1 + S-M to denote Algorithm 2.Algorithm 2:Computing the Moore-Penrose inverse of the Toeplitz matrix H.
Require: Start with H
∈
C
m
×
n, n
>
m. Let h
i be the ith column of H, for i
=
1
,
…
,
n.
(1) (Improved HSR1 Step) DefineE and L as in Theorem 6. Then, set H
m
=
E
-
L.
(2) Sety
m
,
t
=
H
m
-
1
h
t
∈
C
m for all t
=
1
,
…
,
n and X
m
=
H
∗
H
m
-
1.
(3) (S-M Step) Computey
l
,
t and X
l by using (22) and (23), respectively for all l
=
m
+
1
,
…
,
n.
(4) OutputH
†
=
X
n.Note that proposed Algorithms1 and 2 are not using directly, at any step, the Sherman-Morrison formula to compute the inverse of some matrix. In Section 4.4, we provide more information regarding this situation.


## 4.3. Numerical Experiments

In this subsection we analyze numerical data arising during the computation of the Moore-Penrose inverse of the Toeplitz matrixH by applying aMatlab implementation of Algorithms 1 and 2. In order to test the time efficiency as well as the accuracy of considered methods, we enrich our collection of matrices used in Tables 1 and 2 with larger matrices. In order to complete our numerical study, we also compare the results derived by applying Algorithm  2 from [23] (Algorithm BP). A comparison of Algorithms 1 and 2 and BP is presented in Table 3.Table 3
Comparison of Algorithm1 (HSR1 + S-M), Algorithm 2 (IHSR1 + S-M), and BP.
Method
Sizesn, ℓ, s
CPU time
r
1
r
2
r
3
r
4
HSR1 + S-M
50, 15, 10
0.0054
1.6449e − 14
2.5845e − 12
3.5344e − 13
1.9362e − 14
IHSR1 + S-M
50, 15, 10
0.0035
2.0486e − 14
2.7318e − 12
5.4426e − 13
2.3222e − 14
BP
50, 15, 10
0.0018
2.0011e − 15
2.0385e − 14
3.773e − 15
1.6463e − 14
HSR1 + S-M
50, 20, 10
0.0040
1.8748e − 14
8.8674e − 13
2.2185e − 13
1.5379e − 14
IHSR1 + S-M
50, 20, 10
0.0032
9.4733e − 15
8.9539e − 13
1.9076e − 13
1.4775e − 14
BP
50, 20, 10
0.0027
9.2684e − 16
1.6646e − 14
3.2679e − 15
5.8648e − 15
HSR1 + S-M
50, 20, 500
0.0053
1.2862e − 14
9.6327e − 13
1.2288e − 13
1.4447e − 14
IHSR1 + S-M
50, 20, 500
0.0044
8.6971e − 15
1.5068e − 12
2.2331e − 13
1.5216e − 14
BP
50, 20, 500
0.0022
6.8218e − 16
2.1407e − 14
3.2663e − 15
6.0107e − 15
HSR1 + S-M
250, 15, 500
0.0482
2.5535e − 13
3.7054e − 11
4.7471e − 12
2.5608e − 13
IHSR1 + S-M
250, 15, 500
0.0195
2.209e − 13
3.7065e − 11
5.5401e − 12
2.4239e − 13
BP
250, 15, 500
0.0032
3.4724e − 15
2.3693e − 13
1.5467e − 14
7.4894e − 14
HSR1 + S-M
400, 20, 500
0.3234
4.5138e − 13
6.0769e − 10
1.4669e − 11
5.5243e − 13
IHSR1 + S-M
400, 20, 500
0.1466
5.881e − 13
6.0878e − 10
1.5866e − 11
5.3143e − 13
BP
400, 20, 500
0.0116
6.5483e − 15
6.052e − 13
3.0382e − 14
1.3117e − 13
HSR1 + S-M
1200, 20, 500
11.0062
3.0653e − 12
1.7107e − 08
1.4049e − 10
3.2612e − 12
IHSR1 + S-M
1200, 20, 500
2.3018
3.5064e − 12
1.711e − 08
1.3136e − 10
3.1795e − 12
BP
1200, 20, 500
0.1049
1.2805e − 14
3.179e − 12
8.9833e − 14
2.3321e − 13
HSR1 + S-M
1500, 20, 500
21.7539
4.3812e − 12
3.2457e − 08
2.2848e − 10
4.2712e − 12
IHSR1 + S-M
1500, 20, 500
3.7358
4.2345e − 12
3.2456e − 08
2.0156e − 10
4.2919e − 12
BP
1500, 20, 500
0.1886
1.3786e − 14
4.0422e − 12
1.1081e − 13
3.0133e − 13
HSR1 + S-M
2000, 50, 500
74.5230
5.5113e − 12
4.6497e − 08
5.1185e − 10
6.7183e − 12
IHSR1 + S-M
2000, 50, 500
30.8622
5.3682e − 12
4.6504e − 08
5.4329e − 10
6.6102e − 12
BP
2000, 50, 500
0.4038
1.8118e − 14
6.5462e − 12
1.408e − 13
3.6715e − 13From Table3, it is observable that the level of improvement in CPU time that can be achieved by using Algorithm 2 instead of Algorithm 1 is significant. Also, it is clear that the accuracy of both algorithms is almost the same. Clearly, the BP method requires minimal CPU time.


## 4.4. Rounding Error Analysis

The starting motivation of this subsection is the following basic problem: approximation error grows with repeated use of the S-M formula. As a consequence, computations which involve the S-M rule become unstable. For this purpose, it is interesting to investigate the error curves corresponding to the four Penrose equations during the iterations included in step  (3) of Algorithm2. Nevertheless, it is important to note that the proposed Algorithm 2 does not use in each step the S-M rule. Of course, Theorem 5 is in the heart of Algorithm 2 but also it is clear that the poor stability of the S-M formula has no serious influence on the stability of proposed algorithms.So, in this subsection, we shall present error estimations regarding Algorithm2. That is, we record and investigate the residual matrix norms r
1
,
r
2
,
r
3, and r
4, corresponding to the four Penrose equations during the execution of step  (3) of Algorithm 2. The results are presented in Figures 1, 2, and 3 for the cases: (a) n
=
1200, l
=
20, and s
=
500, (b) n
=
1500, l
=
20, and s
=
500, and (c) n
=
2000, l
=
20
,
a
n
d
s
=
500, respectively. The choice of these matrices from Table 3 is random.Figure 1
Errors recorded forn
=
1200, l
=
20, and m
=
500 while the loop of step (3) is iterating.Figure 2
Errors recorded forn
=
1500, l
=
20, and m
=
500 while the loop of step (3) is iterating.Figure 3
Errors recorded forn
=
2000, l
=
50, and m
=
500 while the loop of step (3) is iterating.The purple line in each of Figures1–3 corresponds to the values of the matrix norm (
X
l
H
)
∗
-
X
l
H
2, l
=
m
+
1
,
…
,
n. Since X
l
=
H
∗
H
l
-
1 and the matrix H
l is symmetric, the matrix X
l
H is symmetric. For this reason, the graph corresponding to these values is almost a constant line.The stable convergence is observable from each of these figures. Also, it is possible to notice very fast convergence in the terminal phase of the convergence as well as a relatively slower convergence in the middle of the recurrent process. This fact is understandable if one takes into account that it is necessary to take into consideration all the columns of the input matrix in order to get the complete information on the pseudoinverse.


## 5. Application in Image Restoration

Several image restoration examples are presented in this subsection. Experiments are done using Matlab programming package on an Intel(R) Core(TM) i5-2540M CPU @ 2.6 GHz 64-bit system with 8 GB RAM memory.Numerical results corresponding to the Moore-Penrose inverse are derived using Algorithm2. Results are derived using n
=
1200 and s
=
100 and different values of l. Suppose that the matrix F
∈
R
r
×
m corresponds to the original image and G
∈
R
r
×
m is the matrix corresponding to the degraded image.The model of the nonuniform blurring is the same as in [23] and assumes that the blurring of columns in the image is independent with respect to the blurring of its rows. Therefore, the relations between the original and the blurred image are expressed by the matrix equation(37)
G
=
H
C
F
H
R
T
,
G
∈
R
r
×
m
,
H
C
∈
R
r
×
n
,
F
∈
R
n
×
t
,
H
R
∈
R
m
×
t
,where n
=
r
+
l
c
-
1, t
=
m
+
l
r
-
1, l
c is length of the vertical blurring, and l
r is length of the horizontal blurring (in pixels). Following the approach used in [23], the Moore-Penrose inverse approach is used to restore the blurred image G, which produces the next approximation of F:(38)
F
~
=
H
C
†
G
H
R
†
T
.The algorithm for computing the Moore-Penrose inverse of the blurring matrix used in [24, 25] is based on the method for computing the Moore-Penrose inverse of a full rank matrix, introduced in [26, 27]. In our case, we use Algorithm 2.The first approach is based on the usage of the pure SR1 update algorithm, proposed in [15]. The second approach starts with the SR1 update algorithm and, after its completion, continues with the BP method. The third approach is a hybrid combination which comprises the SR1 updates in the first step and the S-M formula in the second phase (Algorithm 1). The fourth approach is based on the improvement of the SR1 update step method of Algorithm 1, which leads to Algorithm 2 with the best performances. The improvement is achieved using the specific structure of the underlying Toeplitz matrix.The most popular algorithm for image restoration based on the matrix pseudoinverse is developed using the singular value decomposition (SVD) [28]. But, the SVD is sensitive to singular values very close to zero. For this purpose, it is common approach to leave out small singular values, corresponding to high-frequency components. The resulting method is known as the truncated SVD (or TSVD shortly) [28]. On the other hand, the methods proposed in the present paper are based on different approach and do not require any kind of the regularization. In order to verify the applicability of the proposed method in the image restoration, we present a comparison of results generated by the TSVD and the rank-one updates.Parameters of Barbara image degradation arel
r
=
39, l
c
=
25, s
=
l
r
/
2, and s
=
l
c
/
2 and salt pepper noise with the noise density of 0.05 is assumed. The original and damaged image are presented in Figure 4.Figure 4
(a) Image of Barbara. (b) Blurred noisy image of Barbara.
(a)
(b)Comparison with the state-of-the-art methods in image processing is presented. Many of them are constrained least-squares filter (CLS); Wiener filter (WF); Lucy-Richardson algorithm (LR); and truncated singular value decomposition (TSVD) [28]. Improvements in signal-to-noise ratio (ISNR) and Dice Coefficient (DC) are considered.Results corresponding to ISNR values are presented in Figures5 and 6.Figure 5
ISNR versusl
c for the removal of Gaussian white noise with zero mean and variance 0.05 and Gaussian blur (l
r
=
35, s
=
l
r
/
2, and s
=
l
c
/
2) for Barbara image.Figure 6
ISNR versusl
c for the removal of Gaussian white noise with zero mean and variance 0.05 and Gaussian blur (l
r
=
35, and s
=
100) for Barbara image.Values corresponding to Dice Coefficient (DC) [29, 30] are presented in Figure 7.Figure 7
DC versusl
c for the removal of Gaussian white noise with zero mean and variance 0.05 and Gaussian blur (l
r
=
35, and s
=
100) for Barbara image.Restored images are presented in Figures8 and 9.Figure 8
(a) Moore-Penrose restored image. (b) Lucy-Richardson restored image.
(a)
(b)Figure 9
(a) Wiener filter restored image, (b) constrained least-squares filter restored image, and (c) truncated singular value decomposition restored image.
(a)
(b)
(c)


## 6. Conclusion

The present paper investigates the application of the SR1 updates procedure and S-M formula in computation of the Moore-Penrose inverse of specific Toeplitz matrices that appear in the image restoration process.The first approach is based on the usage of the pure SR1 update algorithm, proposed in [15]. The second approach starts with the SR1 update algorithm and, after its completion, continues with the BP method. The third approach is a hybrid combination which comprises the SR1 updates in the first step and the S-M formula in the second phase (Algorithm 1). The fourth approach is based on the improvement of the SR1 update step method of Algorithm 1, which leads to Algorithm 2 with the best performances. The improvement is achieved using the specific structure of the underlying Toeplitz matrix.An application of Algorithm2 in image restoration is considered. In this case, a comparison with the modified BP method from [23] is presented.



---

*Source: 9065438-2016-08-23.xml*