УДК 519.

65
Yu. Hakopian
DOI: https://doi.org/10.18523/2617-70802201911-23

COMPUTING THE MOORE-PENROSE INVERSE FOR

BIDIAGONAL MATRICES
The Moore-Penrose inverse is the most popular type of matrix generalized inverses which has many
applications both in matrix theory and numerical linear algebra. It is well known that the Moore-Penrose
inverse can be found via singular value decomposition. In this regard, there is the most effective algo-
rithm which consists of two stages. In the first stage, through the use of the Householder reflections, an
initial matrix is reduced to the upper bidiagonal form (the Golub-Kahan bidiagonalization algorithm).
The second stage is known in scientific literature as the Golub-Reinsch algorithm. This is an itera-
tive procedure which with the help of the Givens rotations generates a sequence of bidiagonal matrices
converging to a diagonal form. This allows to obtain an iterative approximation to the singular value
decomposition of the bidiagonal matrix.
The principal intention of the present paper is to develop a method which can be considered as an
alternative to the Golub-Reinsch iterative algorithm. Realizing the approach proposed in the study, the
following two main results have been achieved. First, we obtain explicit expressions for the entries of
the Moore-Penrose inverse of bidigonal matrices. Secondly, based on the closed form formulas, we get
a finite recursive numerical algorithm of optimal computational complexity. Thus, we can compute the
Moore-Penrose inverse of bidiagonal matrices without using the singular value decomposition.
Keywords: Moor-Penrose inverse, bidiagonal matrix, inversion formula, finite recursive algorithm.

Introduction The most effective procedure to compute the

Moore-Penrose inverse involves two main stages
As is known, for a real m × n matrix A the [4].
Moore-Penrose inverse A+ is the unique matrix Stage 1. Matrix reduction to the bidiagonal
that satisfies the following four properties [1]: form.
At this stage an m × n matrix, where m ≥ n,
AA+ A = A , A+ AA+ = A+ , by means of the Householder reflections is trans-
formed to upper bidiagonal form
(A+ A)T = A+ A , (AA+ )T = AA+ .  
a11 a12 0 ... 0
If A is a square nonsingular matrix, then A+ =  0 a22 a23 ... 0 
= A−1 . Thus the Moore-Penrose inverse general-  ..

. . . .

izes the ordinary matrix inversion.  . .. .. .. .. 

.
 
There is well-known formula for the Moore-  0 ... 0 an−1 n−1 an−1 n 
 
Penrose inverse which is obtained by the singular  0 ... 0 0 ann  

value decomposition (abbreviated SVD) of the ma-  
trix (see [1; 4], for instance). 0
The singular value decomposition of an m × n (1.3)
matrix A with rank r is its factorization of the form The computational process is known as Golub-
Kahan bidiagonalization [2]. Thereby the problem
A = U ΛV T , (1.1) is reduces to the Moore-Penrose inversion of the
bidiagonal matrix (1.3).
where U is an m × m orthogonal matrix, Λ =
Stage 2. Golub-Reinsch SVD iterative algo-
= diag [σ1 , σ2 , . . . , σr ] is an m × n diagonal matrix,
rithm.
and V is an n×n orthogonal matrix. The diagonal
Once the bidiagonalization of the initial ma-
entries σ1 ≥ σ2 ≥ . . . ≥ σr > 0 of Λ are known as
trix has been achieved, the next task is to zero the
singular values of the matrix A. Having the fac-
superdiagonal entries in the matrix (1.3). With
torization (1.1), the Moore-Penrose inverse can be
this purpose the Golub-Reinsch algorithm is im-
written as
plemented [3]. The algorithm, with the help of the
A+ = V Λ+ U T , (1.2)
Givens rotations generates a sequence of bidiago-
where Λ+ = diag [σ1−1 , σ2−1 , . . . , σr−1 ] is n × m di- nal matrices that converge to a diagonal form. As
agonal matrix. a result, at a certain step of the iterative process we

c Yu. Hakopian, 2019
12 ISSN 2617-7080. Могилянський математичний журнал. 2019. Том 2

get an approximation to the SVD of the bidiagonal Case 1: a11 6= 0, ann 6= 0.

matrix (1.3). Having the SVD, the Moore-Penrose Let zero diagonal entries of the matrix A are
inverse of the matrix is computed (see [1; 4], for ai1 i1 , ai2 i2 , . . . , aip−1 ip−1 , where 1 < i1 < i2 <
instance). < · · · < ip−1 < n and p > 1. We split the matrix
The objective of the present work is to de- into blocks drawing dividing lines after the rows
velop a method which allows to deduce formulas i1 − 1, i2 − 1, . . . , ip−1 − 1 and after the columns
for the entries of the Moore-Penrose inverse of i1 , i2 , . . . , ip−1 . As a result, the matrix (2.1) takes
upper bidiagonal matrices. The obtained closed a block diagonal form. The first and the last diag-
form solution to the Moore-Penrose inversion may onal blocks are rectangular bidiagonal matrices of
be considered as an alternative to sufficiently the sizes (i1 −1)×i1 and (n−ip−1 +1)×(n−ip−1 ),
labour-consuming Golub-Reinsch iterative proce- respectively, while the remaining blocks are square
dure briefly described in the Stage 2 of this sec- lower bidiagonal matrices. As an illustration, a
tion. Moreover, explicit expressions for the entries pattern of the matrix (for n = 10), the partition-
of the Moore-Penrose inverse lead to fairly simple ing procedure and resulting block diagonal struc-
finite numerical algorithm with optimal volume of ture are shown in Fig. 2.1 (stars represent nonzero
computational expenditures. entries).

Partition of a bidiagonal matrix into blocks

? ?
? ?
Let us consider a real n × n upper bidiagonal
? ?
matrix 0 ?
?
 
a11 a12 0
 a22 a23 0  ? ?
? ?
 
A=
 . .. . .. .

  0 ?
 0 an−1 n−1 an−1 n  ? ?
ann ?
(2.1)
Note that it suffices to consider square upper bidi- ? ?
agonal matrices since for rectangular upper bidi- ? ?
agonal matrices the problem can be easily reduced ? ?
to our case. Indeed, if m > n then according to ?
(1.3) we have the block structure ?
? ?
A ? ?
,
0 ?
where A is a square upper bidiagonal matrix of the ? ?
?
form (2.1). It can be seen that in this case
+ Figure 2.1. Partition of the matrix (case 1).
A
= A+ 0T .

0 Case 2: a11 = 0, ann 6= 0.
We allocate the first column of the matrix
We assume that the matrix A is singular, i.e. A, as a separate zero block of the size n ×
a11 a22 . . . ann = 0. Next, we assume that × 1. Next, we partition the remaining subma-
ai i+1 6= 0, i = 1, 2, . . . , n − 1. (2.2) trix into diagonal blocks as follows. If there are
other zero diagonal entries of the matrix A, say
Otherwise, if some of superdiagonal entries of the ai1 i1 , ai2 i2 , . . . , aip−1 ip−1 , where 1 < i1 < i2 <
matrix A are equal to zero, the problem of comput- < · · · < ip−1 < n and p > 1, then the subma-
ing the Moore-Penrose inverse is decomposed into trix is subdivided according to the rule drscribed
several similar problems for bidiagonal matrices of in the Case 1. As an illustration, see a pattern
lower order. of the matrix given in Fig. 2.2 . The last diagonal
To compute the Moore-Penrose inverse of the block of the submatrix is rectangular bidiagonal
matrix A, we apply a special partition of this ma- matrix of the size (n − ip−1 + 1) × (n − ip−1 ); the
trix into blocks. The partitioning procedure uses remaining diagonal blocks are square lower bidiag-
the arrangement of zeros on the main diagonal onal matrices. If there are no other zero diagonal
of the matrix. We distinguish the following four entries, except the first one, then the submatrix is
cases. not subdivided.
Yu. Hakopian. Computing the Moore-Penrose inverse for bidiagonal matrices 13

0 ? ? ?
? ? ? ?
? ? ? ?
? ? ?
0 ? ? ?
0 ? ?
? ? ?
0 ? ? ?
? ? ? ?
? 0

? Figure 2.3. Partition of the matrix (case 3).

? ?
? ?
? ? Case 4: a11 = 0, ann = 0.
? The allocation of the first column and the last
0
? row of the matrix A gives us three zero blocks of
? ? the sizes (n − 1) × 1, 1 × (n − 1) and 1 × 1 (see
? Fig. 2.4). Then we partition the remaining subma-
? ? trix. If there are other zero diagonal entries of the
? matrix A, say ai1 i1 , ai2 i2 , . . . , aip−1 ip−1 , where 1 <
< i1 < i2 < · · · < ip−1 < n and p > 1, then the
submatrix is subdivided by the rule drscribed in
the Case 1. The diagonal blocks of this subdivi-
Figure 2.2. Partition of the matrix (case 2). sion are square lower bidiagonal matrices. If there
are no other zero diagonal entries except the first
and last, then the submatrix is not subdivided.
Case 3: a11 6= 0, ann = 0.
First, we allocate the last row of the matrix
A, as a separate zero block of the size 1 × n. 0 ?
Next, we partition the remaining submatrix into ? ?
diagonal blocks using the same idea. If there are 0 ?
other zero diagonal entries of the matrix A, say ? ?
? ?
ai1 i1 , ai2 i2 , . . . , aip−1 ip−1 , where 1 < i1 < i2 <
0 ?
< · · · < ip−1 < n and p > 1, then the submatrix
0 ?
is subdivided by the rule drscribed in the Case 1 ? ?
(see a pattern of the matrix given in Fig. 2.3). The ? ?
first diagonal block of the submatrix is rectangu- 0
lar bidiagonal matrix of the size (i1 − 1) × i1 ; the
remaining diagonal blocks are square lower bidiag-
?
onal matrices. If there are no other zero diagonal
? ?
entries, except the last one, then the submatrix is ?
not subdivided. ? ?
0 ? ?
?
?
? ? ? ?
? ? ? ?
? ? 0 0
0 ?
? ?
0 ? Figure 2.4. Partition of the matrix (case 4).
0 ?
? ? Thus we have four principal cases of block par-
? ? titioning the initial upper bidiagonal matrix A,
0 schematically presented in Fig. 2.5.
14 ISSN 2617-7080. Могилянський математичний журнал. 2019. Том 2

B1
B1+
B2
B2+
A= · A+ =
·
· ·
·
·
Bp Bp+

case 1 case 1

0
B1
B1+
B2
B2+
A= 0 · A+ =
·
· ·
·
·
Bp
Bp+
case 2 case 2

B1

B2
B1+
A= ·
·
·
B2+
Bp A+ = 0
·
0 ·
·
case 3 Bp+

B1 case 3

B2 0 0

A= 0 · B1+
·
·
B2+
Bp A+ = 0
·
0 0 ·
·
case 4 Bp+

case 4

Figure 2.5. The cases of block partitioning.

Accordingly, the Moore-Penrose inverse also Figure 2.6. The structure of the Moore-Penrose
has a block structure, as shown in Fig. 2.6. inverse.
Yu. Hakopian. Computing the Moore-Penrose inverse for bidiagonal matrices 15
+
Summarizing the previous reasoning, we con- Proposition 1. The entries of the matrix B =
clude that our task is to find the Moore-Penrose = [zij ]m×m are as follows: for the indices i =
inverses for blocks of the following three types: = 1, 2, . . . , m we have
type 1 : bidiagonal block of a size m × m; i−1
1 Y
type 2 : bidiagonal block of a size m − 1 × m; zij = (−1)i+j rs , j = 1, 2, . . . , i − 1;
di s=j
type 3 : bidiagonal block of a size m × m − 1.
In Fig. 2.7 we schematically give the types of di- zii = d1i ;
agonal blocks (the mark ? stands for a nonzero zij = 0, j = i + 1, i + 2, . . . , m,
entry). (2.4)
where
?
? ? bs
rs ≡ , s = 1, 2, . . . , m − 1. (2.5)
? ? ds
.. .. Based on the formulas (2.4) we can write the
. .
following simple procedure to calculate the entries
? ? of the matrix B + .
type 1 Algorithm (B ⇒ B+ )/type 1
1. Compute the quantities rs defined in (2.5).
? ? 2. Compute the lower triangular part of the
? ? matrix B + . For indices i = 1, 2, . . . , m:
.. .. 1
. . zii = ; zij = −rj zi j+1 , j = i − 1, i − 2, . . . , 1.
? ? di
End algorithm
type 2
It can be readily seen that Algorithm
(B ⇒ B+ )/type 1 requires
?
1 2
? ? A(1)
ops = m + O(m) (2.6)
.. 2
? . arithmetical operations.
.. Next, we will focus our attention on computing
. ?
? the Moore-Penrose inverse for the blocks of type 2
and 3.
type 3
A way of computing the Moore-Penrose
Figure 2.7. The types of diagonal blockfs. inverse
It is necessary to pay attention to the follow- To solve the problem, in this section we outline
ing circumstance. As follows from the process of an approach based upon the well-known equality
partitioning the initial matrix (2.1) into blocks, in
each of the Cases 1-4 we have at most two rect- A+ = lim (AT A + εI)−1 AT , (3.1)
ε→+0
angular blocks (of size m − 1 × m or m × m − 1).
The remaining blocks are square lower bidiagonal where I is the identity matrix, which holds true
matrices. As an illustration, see Fig. 2.5. for any real matrix (see [1; 4], for instance). Here
we present the main ideas to compute the Moore-
Computing the Moore-Penrose inverse for a
Penrose inverse for a block of the type 2. For a
block of the type 1 is not difficult. Consider a
block of the type 3, as will be seen below, the
square matrix
problem is reduced to the case under considera-
  tion. Note that, as in the previous case, it is con-
d1
 b1 venient to introduce new notation for the entries
d2 
B=

.. .. ,

(2.3) of the block.
 . .  Let us have an m − 1 × m bidiagonal matrix
bm−1 dm 
d1 b1

 d2 b2 
where di 6= 0, i = 1, 2, . . . , m and bi 6= 0, i = 
.. ..

B= ,
 
= 1, 2, . . . , m−1 (we choose new notation for block  . . 
entries). Since the matrix (2.3) is nonsingular then  dm−2 bm−2 
B + = B −1 (see [1], for instance). This inverse can dm−1 bm−1
be easily found. (3.2)
16 ISSN 2617-7080. Могилянський математичний журнал. 2019. Том 2

where bi , di 6= 0, i = 1, 2, . . . , m − 1 and m ≥ 2. is tridiagonal matrix of the following structure:

The matrix
L(ε) ≡ B T B + εI (3.3)
 2 
d1 + ε b 1 d1
 b1 d1 d22 + b21 + ε b2 d2 0 
 
L(ε) = 
 . .. . .. . .. .

(3.4)
 
2 2
 0 bm−2 dm−2 dm−1 + bm−2 + ε bm−1 dm−1 
bm−1 dm−1 b2m−1 + ε
To invert this matrix, let us apply advantage of 6. The entries of the lower triangular part of
the algorithm developed in [5]. Consider a nonsin- the matrix C −1 are found:
gular symmetric tridiagonal matrix
xij = xji , i = j + 1, j + 2, . . . , m ;
 
c11 c12
 c21 c22 j = 1, 2, . . . , m − 1 . (3.11)
 c23 0 

C= .. .. .. End procedure
.
 
. . .
  The proposed way to obtain the Moore-Penrose
 0 cm−1 m−2 cm−1 m−1 cm−1 m  inverse B + is as follows. In consistence with equal-
cm m−1 cm m ity (3.1) and notation (3.3), we have
(3.5)
−1
We assume that m ≥ 2. Referring to [5], the ma- B + = lim L(ε) B T . (3.12)
ε→+0
trix C −1 = [xij ]m×m can be obtained by the fol-
lowing computational procedure. −1
Finding first the inverse matrix L(ε) , the entries
−1
Procedure 3d/inv (C ⇒ C −1 ) of the matrix L(ε) B T are calculated and a char-
1. Compute the quantities fi (i = 2, 3, . . . , m), acter of their dependence on the parameter ε is
gi (i = 2, 3, . . . , m − 1) revealed. Then, according to the equality (3.12),
and hi (i = 1, 2, . . . , m − 1): passing to the limit when ε → +0, we arrive to a
closed form expressions for the entries of the ma-
cii ci i+1 cii +
fi = , gi = , hi = . (3.6) trix B .
ci i−1 ci i−1 ci i+1
Note: if m = 2, then the quantities The Moore-Penrose inverse of rectangular
gi are not introduced. blocks
2. Compute recursively the quantities µi (i =
= 1, 2, . . . , m): Let us consider as the matrix C from (3.5) the
tridiagonal matrix L(ε) obtained in (3.4). Com-
µm = 1 , µm−1 = −fm , paring the records of these matrices, we have
µi = −fi+1 µi+1 − gi+1 µi+2 , (3.7) cii = d2i + b2i−1 + ε , i = 1, 2, . . . , m (4.1)
i = m − 2, m − 3, . . . , 1.
(in order to unify records of formulas, we set b0 =
3. Compute recursively the quantities νi (i = = 0) and
= 1, 2, . . . , m):
ci i+1 = bi di , i = 1, 2, . . . , m − 1 ;
ν1 = 1 , ν2 = −h1 , ci i−1 = bi−1 di−1 , i = 2, 3, . . . , m . (4.2)
1 In accordance with our plan, let us carry out
νi = −hi−1 νi−1 − νi−2 , i = 3, 4, . . . , m.
gi−1 a more detailed analysis of the quantities succes-
(3.8) sively computed in the procedure 3d/inv from
4. Compute the quantity Section 3.
Consider first the quantities fi , gi and hi which
t = (c11 µ1 + c12 µ2 )−1 . (3.9)
were introduced in (3.6). Using the expressions
Note: since C is nonsingular matrix, (4.1) and (4.2), we get
then c11 µ1 + c12 µ2 6= 0 [5]. ◦
5. The entries of the upper triangular part of fi =f i +αi ε , i = 2, 3, . . . , m,
the matrix C −1 are computed:
where
xij = µj νi t, i = 1, 2, . . . , j ; j = 1, 2, . . . , m . ◦ d2i + b2i−1 1
(3.10) f i= , αi = ; (4.3)
bi−1 di−1 bi−1 di−1
Yu. Hakopian. Computing the Moore-Penrose inverse for bidiagonal matrices 17
bi di ◦
gi = , i = 2, 3, . . . , m − 1; (4.4) Lemma 3. The quantities µi can be written in
bi−1 di−1 the form
◦
hi =hi +βi ε , i = 1, 2, . . . , m − 1, m−1
◦ Y
where µi = (−1)m−i rs , i = 1, 2, . . . , m. (4.10)
d2i + b2i−1
◦ 1 s=i
hi = , βi = . (4.5)
bi di bi di ◦
Next, go to the quantities µi and νi recursively Proof. Firstly, the value µm = 1 conforms to
defined in (3.7) and (3.8), respectively. the record (4.10). Then, in accordance with (4.3)
Lemma 2. The quantities µi are represented as and (4.7),
◦ ◦ ◦ ◦ bm−1
µm =µm +γm ε , µm−1 =µm−1 +γm−1 ε , µm−1 = − f m = − = −rm−1 .
dm−1
◦
µi =µi +γi ε + O(ε2 ) , 1 ≤ i ≤ m − 2,
Further reasoning is carried out by induction. Us-
(4.6) ing the expressions (4.3) and (4.4), proceeding
◦
where the quantities µi and γi satisfy the following from (4.7) we obtain
recurrence relations:
m−1
◦ ◦ ◦ ◦ d2i+1 + b2i Y
µm = 1 , µm−1 = − f m , µi = − (−1)m−i−1 rs
bi di s=i+1
◦ ◦ ◦ ◦ (4.7) m−1
µi = − f i+1 µi+1 −gi+1 µi+2 , bi+1 di+1 Y
i = m − 2, m − 3, . . . , 1 − (−1)m−i−2 rs
bi di s=i+2
and
m−1
γm = 0 , γm−1 = −αm , d2i+1 + b2i

m−i
Y bi+1 di+1
◦ ◦
= (−1) rs ri+1 −
bi di bi di
γi = − f i+1 γi+1 − gi+1 γi+2 − αi+1 µi+1 , s=i+2
i = m − 2, m − 3, . . . , 1. m−1
Y
(4.8) = (−1)m−i ri ,
◦
Proof. Since µm = 1 then in (4.6) we set µm = s=i
= 1, γm = 0. Further, µm−1 = −fm (see (3.7)). which completes the proof. 2
According to the expressions (4.3) we have fm = The next assertion is a simple consequence of
◦
=f m +αm ε. Therefore in the representation (4.6) the formula (4.10).
◦ ◦ Corollary 4. The following relation holds:
we set µm−1 = − f m , γm−1 = −αm .
◦ ◦
For the indices in the range 1 ≤ i ≤ m − 2, µi = −ri µi+1 , i = 1, 2, . . . , m − 1 . (4.11)
required representations can be readily derived by
induction from the relations (3.7) using expressions A representation similar to (4.6) takes place
(4.3). Indeed, having done simple transformations also for the quantities νi .
as follows Lemma 5. The quantities νi are represented as
µi = −fi+1 µi+1 − gi+1 µi+2 ◦ ◦
◦
ν1 =ν 1 +δ1 ε , ν2 =ν 2 +δ2 ε ,
◦
= −(f i+1 +αi+1 ε)(µi+1 +γi+1 ε + O(ε2 )) ◦
(4.12)
◦ νi =ν i +δi ε + O(ε2 ) , 3 ≤ i ≤ m,
−gi+1 (µi+2 +γi+2 ε + O(ε2 ))
◦ ◦ ◦ ◦
= (− f i+1 µi+1 −gi+1 µi+2 ) where the quantities ν i and δi satisfy the following
◦ ◦ recurrence relations:
+(− f i+1 γi+1 − gi+1 γi+2 − αi+1 µi+1 )ε
◦ ◦ ◦
+O(ε2 ) , ν 1 = 1 , ν 2 = − h1 ,
◦ ◦ ◦ 1 ◦
we get (4.6) as well as recurrence relations (4.7) ν i = − hi−1 ν i−1 − ν i−2 , i = 3, 4, . . . , m
and (4.8). 2 gi−1
◦ (4.13)
The quantities µi computed by the recursion and
(4.7) can be represented in closed form.
Let us introduce the following notation: δ1 = 0 , δ2 = −β1 ,
◦ 1 ◦
bs δi = − hi−1 δi−1 − δi−2 − βi−1 ν i−1 ,
rs ≡ , s = 1, 2, . . . , m − 1 . (4.9) gi−1
ds i = 3, 4, . . . , m.
Additionally, we set r0 = rm = 1. (4.14)
18 ISSN 2617-7080. Могилянський математичний журнал. 2019. Том 2
◦
Proof. Since ν1 = 1 then in (4.12) we set ν 1 = 1, Corollary 7. The following relation holds:
δ1 = 0. Further, ν2 = −h1 (see (3.8)). Accord- ◦ 1 ◦
◦ ν i+1 = − νi , i = 1, 2, . . . , m − 1 . (4.16)
ing to the expressions (4.5) we have h1 =h1 +β1 ε. ri
◦
Therefore in the representation (4.12) we set ν 2 = Our next task is to derive an expression for the
◦
= − h1 , δ2 = −β1 . quantity t given in (3.9), depending on the param-
For the indices in the range 3 ≤ i ≤ m, re- eter ε. Since c11 = d21 +ε, c12 = b1 d1 (see (4.1) and
quired representations can be readily derived by (4.2)) then taking into account the representations
induction from the relations (3.8) using expressions (4.6) for the quantities µi we get
(4.5). Indeed, having done simple transformations ◦
t = ((d21 + ε)(µ1 +γ1 ε + O(ε2 ))
as follows ◦
+b1 d1 (µ2 +γ2 ε + O(ε2 )))−1
1 ◦ ◦
νi = −hi−1 νi−1 − νi−2 = (d21 (µ1 +r1 µ2 )
gi−1 ◦
◦ ◦ +(µ1 +d1 (γ1 d1 + γ2 b1 ))ε + O(ε2 ))−1 .
= −(hi−1 +βi−1 ε)(ν i−1 +δi−1 ε + O(ε2 ))
1 ◦ ◦ ◦
− gi−1 (ν i−2 +δi−2 ε + O(ε2 )) By virtue of the relation (4.11), µ1 +r1 µ2 = 0.
◦ ◦ 1 ◦ Thus
= (− hi−1 ν i−1 − ν i−2 ) 1
gi−1 t= ◦ . (4.17)
◦ ◦
1
+(− hi−1 δi−1 − gi−1 δi−2 − βi−1 ν i−1 )ε (µ1 +d1 (γ1 d1 + γ2 b1 ))ε + O(ε2 )
+O(ε2 ) , Having the representations for the quantities
µi , νi and t, by formulas (3.10) and (3.11) we get
we get (4.12) as well as recurrence relations (4.13)
the entries of the inverse matrix
and (4.14). 2
−1
We can write closed form expressions for the L(ε) = [xij ]m×m .
◦
quantities ν i as well. Further, let us introduce the matrix
◦
Lemma 6. The quantities ν i can be written in −1
the form Y (ε) ≡ L(ε) BT . (4.18)

i−1 According to the equality (3.1) and notation (3.3),

◦ i+1
Y 1 B + = lim Y (ε). If
ν i = (−1) , i = 1, 2, . . . , m. (4.15) ε→+0
r
s=1 s
+
B = [zij ]m×m−1 , Y (ε) = [yij (ε)]m×m−1
◦
Proof. The value ν 1 = 1 conforms to the record then
(4.15). Then in accordance with (4.5) and (4.13), zij = lim yij (ε),
ε→+0
◦ ◦ d 1
ν 2 = − h1 = − 1 = − . i = 1, 2, . . . , m , j = 1, 2, . . . , m − 1. (4.19)
b1 r1
As follows from (4.18), the entries of the matrix
Further reasoning is carried out by induction. Tak- Y (ε) are calculated by the rule
ing into account the expressions (4.4), (4.5) and
yij (ε) = xij dj + xi j+1 bj . (4.20)
using (4.13) we get
Subject to the formulas (3.10) and (3.11), for a
◦ ◦ ◦ 1 ◦
νi = − hi−1 ν i−1 − ν i−2 fixed index j in the range 1 ≤ j ≤ m − 1 we con-
gi−1 sider separately two cases: i = 1, 2, . . . , j and i =
d2i−1 + b2i−2 Y 1 i−2 = j + 1, j + 2, . . . , m.
= − (−1)i • Indices i = 1, 2, . . . , j.
bi−1 di−1 rs
bi−2 di−2 Qs=1
i−3 1
Taking the expression (3.10) for the entries xij ,
i−1
− bi−1 di−1 (−1) s=1 rs from (4.20) we can write
2 2
i+1 di−1 + b i−2 1 bi−2 di−2
= (−1) − yij (ε) = tνi (µj dj + µj+1 bj ). (4.21)
Qi−3 1 bi−1 di−1 ri−2 bi−1 di−1
· s=1 rs Then, using the representations (4.6) of the quan-
1
i−3
1 Y 1 Y 1 tities µi , we have
i−1
= (−1)i+1 = (−1)i+1 , ◦
ri−2 ri−1 s=1 rs r
s=1 i µj dj + µj+1 bj = (µj +γj ε + O(ε2 ))dj
◦
+(µj+1 +γj+1 ε + O(ε2 ))bj
which completes the proof. 2
◦ ◦
The next assertion is a simple consequence of = (µj dj + µj+1 bj )
the formula (4.15). +(γj dj + γj+1 bj )ε + O(ε2 ).
Yu. Hakopian. Computing the Moore-Penrose inverse for bidiagonal matrices 19
As follows from the relation (4.11), With regard of the notation (4.23) we arrive at the
◦ ◦ ◦ ◦
equality
µj dj + µj+1 bj = dj (µj +rj µj+1 ) = 0.
dj+1 1 ◦
uj = − uj+1 − µj+1 . (4.27)
Thus bj bj

µj dj + µj+1 bj = (γj dj + γj+1 bj )ε + O(ε2 ). (4.22) Summarizing the above considerations, on the ba-
sis of the obtained equalities (4.26) and (4.27), we
Substituting the expression (4.22) as well as the can state that the quantities uj satisfy the follow-
representations (4.12) and (4.17) of the quantities ing relations:
νi and t, respectively, into the right hand side of
1
the equality (4.21) yields um−1 = − ,
bm−1
◦ ◦
ν i (γj dj + γj+1 bj ) + O(ε) dj+1 uj+1 + µj+1 (4.28)
yij (ε) = ◦ . uj = − ,
µ1 +d1 (γ1 d1 + γ2 b1 ) + O(ε) bj
j = m − 2, m − 3, . . . , 1 .
By taking limit in the previous equality, according
to (4.19) we find The quantities uj can be represented in closed
◦
form as well. Namely, the following statement
ν i (γj dj + γj+1 bj ) holds.
zij = ◦ , i = 1, 2, . . . , j .
µ1 +d1 (γ1 d1 + γ2 b1 ) Lemma 8. The quantities uj are written as
 
m−j m−k m−1
!
Further, let us introduce the notation (−1)m−j X  Y 1  Y
uj = rs ,
dj r
s=j s
uj ≡ γj dj + γj+1 bj , j = 1, 2, . . . , m − 1 . (4.23) k=1 s=m−k+1

Then the entries zij can be written as follows: j = 1, 2, . . . , m − 1 . (4.29)

◦ The assertion can be proven by direct substitut-
ν i uj ing the expression (4.29) into the relations (4.28)
zij = , i = 1, 2, . . . , j , (4.24)
q and using the expression (4.10) for the quantities
◦
where µj .
◦
q ≡µ1 +d1 u1 . (4.25) As a direct consequence of the expressions
(4.10) and (4.29) we get the expression for the
Now let us turn to the quantities uj defined in
quantity q defined in (4.25).
(4.23). For the index j = m − 1, using the expres-
Lemma 9. The quantity q is written as
sions (4.8) and (4.3), we have
m m−k
! m−1
!
1 X Y 1 Y
m−1
um−1 = γm−1 dm−1 +γm bm−1 = −αm dm−1 = − . q = (−1) rs .
bm−1 r
s=1 s
k=1 s=m−k+1
(4.26) (4.30)
For the indices j = m − 2, m − 3, . . . , 1, taking the Finally, let us replace the expressions (4.15),
expressions (4.8), (4.3) and (4.4), we get ◦
(4.29) and (4.30) of the quantities ν i , uj and q, re-
◦ ◦ spectively, into (4.24). As a result, we obtain the
γj dj = (− f j+1 γj+1 − gj+1 γj+2 − αj+1 µj+1 )dj following expression for the entries of the upper
d2j+1 + b2j bj+1 dj+1 triangular part of the matrix B + :
= − γj+1 dj − γj+2 dj
bj dj bj dj  
◦ m−j
!
− bj1dj µj+1 dj X m−k Y 1 m−1
Y
i+j
(−1)   rs
d2j+1 k=1 s=j s
r
s=m−k+1
= −γj+1 bj − γj+1 zij = i−1 !,
bj m m−k
Y 1
! m−1
bj+1 dj+1 ◦ Y X Y
1 µ
− bj γj+2 − bj j+1 rs · dj rs
s=1 s
r
dj+1 s=1 k=1 s=m−k+1
= −γj+1 bj − (γj+1 dj+1 + γj+2 bj+1 )
bj i = 1, 2, . . . , j . (4.31)
◦
− b1j µj+1 .
• Indices i = j + 1, j + 2, . . . , m.
Hence Using the expressions (3.10) and (3.11), from
(4.20) we get the equality
dj+1 1 ◦
γj dj +γj+1 bj = − (γj+1 dj+1 +γj+2 bj+1 )− µj+1 . yij (ε) = tµi (νj dj + νj+1 bj ). (4.32)
bj bj
20 ISSN 2617-7080. Могилянський математичний журнал. 2019. Том 2

In accordance with the representations (4.12) we Hence

have
◦ bj−1 1 ◦
νj dj + νj+1 bj = (ν j +δj ε + O(ε2 ))dj δj dj +δj+1 bj = − (δj−1 dj−1 +δj bj−1 )− νj .
◦ dj dj
+(ν j+1 +δj+1 ε + O(ε2 ))bj
◦ ◦
= (ν j dj + ν j+1 bj ) In accordance with notation (4.34) we get the
+(δj dj + δj+1 bj )ε + O(ε2 ). equality
As follows from the relation (4.16),
bj−1 1 ◦
◦ ◦ ◦ ◦ wj = − wj−1 − νj . (4.37)
ν j dj + ν j+1 bj = dj (ν j +rj ν j+1 ) = 0. dj dj

Thus Summing up the above considerations, on the ba-

νj dj + νj+1 bj = (δj dj + δj+1 bj )ε + O(ε2 ). (4.33) sis of the equalities (4.36) and (4.37), we infer that
the quantities wj satisfy the following relations:
Substituting the expression (4.33) as well as the
representations (4.6) and (4.17) of the quantities 1
µi and t, respectively, into the right hand side of w1 = − ,
d1
the equality (4.32) yields ◦
bj−1 wj−1 + ν j
◦
µi (δj dj + δj+1 bj ) + O(ε) wj = − , j = 2, 3, . . . , m − 1 .
dj
yij (ε) = ◦ .
µ1 +d1 (γ1 d1 + γ2 b1 ) + O(ε) (4.38)
As with uj , the quantities wj can be repre-
By taking limit in this equality, when ε → +0, sented in closed form. The following statement
according to (4.19) we find holds.
◦ Lemma 10. The quantities wj are written as
µi (δj dj + δj+1 bj )
zij = ◦ , i = j+1, j+2, . . . , m .
µ1 +d1 (γ1 d1 + γ2 b1 ) j k−1
! j−1
!
(−1)j X Y 1 Y
Similarly to the previous case, we introduce the wj = rs ,
dj s=1
rs
k=1 s=k
notation
wj ≡ δj dj + δj+1 bj , j = 1, 2, . . . , m − 1 . (4.34)
j = 1, 2, . . . , m − 1 . (4.39)
Then the entries zij can be written as follows:
◦ The assertion can be proven by direct substi-
µi wj
zij = , i = j + 1, j + 2, . . . , m . (4.35) tution of the expression (4.39) into the relations
q
(4.38) and by using the expression (4.15) for the
◦
Consider the quantities wj defined in (4.34). quantities ν j .
For the index j = 1, using the expressions (4.14) Finally, let us replace the expressions (4.10),
and (4.5), we have ◦
(4.39) and (4.30) of the quantities µi , wj and q,
1 respectively, into the equality (4.35). Resulting
w1 = δ1 d1 + δ2 b1 = −β1 b1 = − . (4.36)
d1 formula for the entries of the lower triangular part
For the indices j = 2, 3, . . . , m − 1, taking the ex- of the matrix B + is of the following type:
pressions (4.14), (4.4) and (4.5) yields
j j−1
m−1
! k−1
! !
◦ 1 ◦ i+j+1
Y X Y 1 Y
δj+1 bj = (− hj δj − δj−1 − βj ν j )bj (−1) rs rs
gj s=i s=1
rs
k=1 s=k
d2j + b2j−1 bj−1 dj−1 zij = m m−k
! m−1
! ,
= − δj bj − δj−1 bj X Y 1 Y
bj dj bj dj dj rs
◦ rs
− bj1dj ν j bj k=1 s=1 s=m−k+1

b2j−1 bj−1 dj−1

= −δj dj − δj − δj−1
dj dj i = j + 1, j + 2, . . . , m . (4.40)
◦
− d1j ν j
bj−1 Combining the above considerations, i.e. hav-
= −δj dj − (δj bj−1 + δj−1 dj−1 )
dj ing the formulas (4.31) and (4.40), we arrive at the
◦
− d1j ν j . following statement.
Yu. Hakopian. Computing the Moore-Penrose inverse for bidiagonal matrices 21
◦
Theorem 11. Let B be an m − 1 × m bidiago- 3. Compute the quantities νi (see
nal matrix given in (3.2). Then the entries of the (4.13),(4.16)):
Moore-Penrose inverse B + = [zij ]m×m−1 of this
◦ ◦ 1 ◦
matrix are as follows: ν 1 = 1 ; ν i+1 = − ν i , i = 1, 2, . . . , m − 1 .
1) for indices j = 1, 2, . . . , m − 1 and i = ri
= 1, 2, . . . , j: 4. Compute the quantities uj (see (4.28)):

m−j
  ! 1
X m−k Y 1 m−1
Y um−1 = − ;
(−1) i+j   rs bm−1
r
s=j s
k=1 s=m−k+1 ◦
zij = i−1 ! !; dj+1 uj+1 + µj+1
m m−k
Y 1 m−1 uj = − ,
bj
Y X Y
rs · dj rs
s=1
r
s=1 s j = m − 2, m − 3, . . . , 1 .
k=1 s=m−k+1
(4.41)
5. Compute the quantities wj (see (4.38)):
2) for indices j = 1, 2, . . . , m − 1 and i = j +
+ 1, j + 2, . . . , m: 1
w1 = − ;
d1
j j−1
m−1
! k−1
! !
i+j+1
Y X Y 1 Y ◦
(−1) rs rs bj−1 wj−1 + ν j
s=i s=1
rs wj = − ,
k=1 s=k dj
zij = m m−k
! m−1
! ,
X Y 1 Y j = 2, 3, . . . , m − 1 .
dj rs
s=1
rs 6. Compute the quantity q (see (4.25)):
k=1 s=m−k+1
(4.42) ◦
where the quantities rs are defined in (4.9). q =µ1 +d1 u1 .
Below is an example to illustrate Theorem 4.1.
7. Compute the upper triangular part of the
Example 1. Consider m − 1 × m bidiagonal matrix
matrix B + (see (4.24)):
 
1 1 ◦
ν i uj
B=
 .. .. .
 zij = , i = 1, 2, . . . , j ; j = 1, 2, . . . , m − 1 .
. . q
1 1
8. Compute the lower triangular part of the
matrix B + (see (4.35)):
Calculations by the formulas (4.41) and (4.42) give
the following result: ◦
µi wj
zij = ,

j
q
i+j

 (−1) 1 − , i = 1, 2, . . . , j ,
m i = j + 1, j + 2, . . . , m ;

zij = ,
 (−1)i+j+1 j , i = j + 1, j + 2, . . . , m j = 1, 2, . . . , m − 1 .


m End algorithm
j = 1, 2, . . . , m − 1 . Note that the numerical implementation of the
Algorithm (B ⇒ B+ )/type 2 requires
Thus in Theorem 4.1 we give formulas for the
A(2) 2
ops = m + O(m) (4.43)
entries of the Moore-Penrose inverse of a block of
the type 2. In addition, based on the expressions arithmetical operations.
and recurrence relations obtained in this section, Next consider a block of type 3. Let an m ×
we suggest a numerical algorithm to compute the × m − 1 bidiagonal matrix
entries of the matrix B + = [zij ]m×m−1 .  
Algorithm (B ⇒ B+ )/type 2 d1
 b1 d2
1. Compute the quantities rs (see (4.9)):

 
 .. 
 b2 . 
bs B=   (4.44)
rs = , s = 1, 2, . . . , m − 1 ; r0 = rm = 1 . . .. d

ds
 m−2

 
 bm−2 dm−1 
◦ bm−1
2. Compute the quantities µi (see (4.7),(4.11)):
◦ ◦ ◦
be given, where bi , di 6= 0, i = 1, 2, . . . , m − 1 and
µm = 1 ; µi = −ri µi+1 , i = m − 1, m − 2, . . . , 1 . m ≥ 2. The problem of finding the Moore-Penrose
22 ISSN 2617-7080. Могилянський математичний журнал. 2019. Том 2

inverse of this matrix is reduced to the previous arithmetical operations.

case. Indeed, by virtue of the well-known property So based on the above study we can formulate
(see [1], for instance) we have the following statement.
B + = ((B T )+ )T , (4.45) Proposition 13. Let a singular upper bidiago-
nal matrix A given in (2.1), with nonzero super-
where the matrix B T is already a block of type 2 diagonal entries, be represented in the block form,
(compare (3.2) and (4.44)). Therefore as a conse- according to the rule described in the Section 2
quence of Theorem 4.1 we can formulate the fol- (Cases 1-4). These are blocks Bk , k = 1, 2, . . . , p
lowing statement. (see Fig. 2.5, as an illustration). Then depend-
Theorem 12. Let B be an m × m − 1 bidiago- ing on the type of a block the entries of the blocks
nal matrix given in (4.44). Then the entries of the B + in block representation of the matrix A+ (see
k
Moore-Penrose inverse B + = [zij ]m−1×m of this Fig. 2.6, as an illustration) are calculated by the
matrix are as follows: formulas obtained in Proposition 2.1 ((2.4) and
1) for indices i = 1, 2, . . . , m − 1 and j = (2.5)), Theorem 4.1 ((4.41) and (4.42)) or The-
= 1, 2, . . . , i: orem 4.2 ((4.46) and (4.47)).
To compute the Moore-Penrose inverse A+ , we
m−i m−k
! m−1
!
X Y 1 Y
i+j
(−1) rs have developed numerical procedures as well.
s=i s
r
zij = j−1
k=1 s=m−k+1 Proposition 14. The entries of the blocks
!;
Bk+ , k = 1, 2, . . . , p, included in the block
m m−k
! m−1
Y X Y 1 Y
rs · di rs structure of the matrix A+ (see Fig. 2.6) can
rs
s=1 k=1 s=1 s=m−k+1 be calculated by Algorithm (B ⇒ B+ )/type 1,
(4.46) Algorithm (B ⇒ B+ )/type 2 or Algorithm
2) for indices i = 1, 2, . . . , m − 1 and j = i + (B ⇒ B+ )/type 3, with expenditure of arithmeti-
+ 1, i + 2, . . . , m: cal operations estimated in (2.6), (4.43) or (4.48),
correspondingly.
  ! !
m−1
Y X i k−1
Y 1 i−1
Y
i+j+1 
(−1) rs  rs Below we give an example to illustrate the work
rs
s=j k=1 s=1 s=k of the numerical algorithms.
zij = m m−k
! m−1
! , Example 2. Consider a matrix, which is divided
X Y 1 Y
di rs into blocks as follows:
s=1 s
r
k=1 s=m−k+1
(4.47) 2 5
 
where the quantities rs are defined in (4.9).  3 −7 
Finally, with equality (4.45) in mind, we can 0 6
 
 
write the following numerical algorithm to com-

 4 2 

+
pute the entries of the matrix B = [zij ]m−1×m .
 −5 −1 
A=  .
Algorithm (B ⇒ B )/type 3 +
 0 4 

T T + 0 2
1. Use Algorithm (B ⇒ (B ) )/type 2 to
 
3 −4
 
T +
compute the matrix (B ) .
 
−6 3
 
+ T + T
2. Calculate B = ((B ) ) . 8
End algorithm
As an obvious consequence of (4.43), we state
Applying the above algorithms, we get the follow-
that Algorithm (B ⇒ B+ )/type 3 also requires
ing result, which coincides with the computations
A(3) 2
ops = m + O(m) (4.48) done with MATLAB.

0.0796 −0.0206
 
0.1682 0.0082
 0.0721 −0.1393 
0.1667 0.0000 0.0000
 
 
+ −0.3333 0.5000 0.0000
A = .
 
1.6667 −2.5000 −1.0000
 0.2500

 
 0.2006 0.1996 −0.1331 0.0499 
0.0506 −0.0337 −0.1442 0.0541
0.0125 −0.0083 0.0055 0.1229

Conclusion

In the work we obtained both explicit formu-

las and finite numerical algorithm to compute the
Yu. Hakopian. Computing the Moore-Penrose inverse for bidiagonal matrices 23
Moore-Penrose inverse of bidiagonal matrices. We (see (2.6), (4.43) and (4.48)), we can assert that
emphasize the following important feature of the for computing one nonzero entry of the matrix
numerical algorithm. Proceeding from the struc- A+ asymptotically one arithmetical operation is
ture of the blocks Bk , k = 1, 2, . . . , p, in the expended. Thereby the proposed computational
block representation of the matrix A+ (namely, the method can be considered as optimal. What is
presence of zeros located at predetermined places) more, we point out another important property of
and the estimations of the number of arithmeti- the computational algorithm. The blocks Bk+ are
cal operations required to compute each block Bk+ computed independently of each other.

References

1. A. Ben-Israel and T. N. E. Greville, Generalized In- position and Least Squares Solutions”. Numer. Math.
verses. Theory and Applications, 2nd ed. (Springer, 14, 403–420, (1970).
New York, 2003). 4. G. H. Golub and Ch. F. van Loan, Matrix Compu-
2. G. H. Golub and W. Kahan., “Calculating the Sin- tations, 3rd ed. (The John Hopkins University Press,
gular Values and Pseudo-Inverse of a Matrix”. SIAM 1996.)
J. Num. Anal., 2, 205–224. (1965) 5. J. W. Lewis, “Invertion of tridiagonal matrices”. Numer.
3. G. H. Golub and C. Reinsch, “Singular Value Decom- Math., 38, 333–345 (1982).

Акопян Ю. Р.

ОБЧИСЛЕННЯ ОБЕРНЕНОГО ВIДОБРАЖЕННЯ

МУРА–ПЕНРОУЗА ДЛЯ ДВУДIАГОНАЛЬНИХ
МАТРИЦЬ
Обернене вiдображення Мура–Пенроуза є найбiльш поширеним вiдображенням, що використо-
вується для пошуку оберненої матрицi. Це вiдображення має численнi застосування як у теорiї
матриць, так i в обчислювальнiй лiнiйнiй алгебрi. Вiдомо, що обернена матриця Мура–Пенроуза
може бути отримана через сингулярний розклад. Найефективнiший з iснуючих алгоритмiв скла-
дається з двох крокiв. На першому кроцi, використовуючи вiдображення Хаусхолдера, початкова
матриця зводиться до верхнього двудiагонального вигляду (алгоритм Голуба–Кахана). Другий
крок вiдомий у науковiй лiтературi як алгоритм Голуба–Райнша. Ця iтерацiйна процедура за
допомогою методу Гiвенса генерує послiдовнiсть двудiагональних матриць, яка збiгається до дiа-
гонального вигляду. В такий спосiб отримується iтерацiйне наближення до сингулярного розкладу
двудiагональної матрицi. Головною метою цiєї статтi є розробка методу, який можна розглядати
як альтернативну замiну алгоритму Голуба–Райнша. За допомогою реалiзацiї запропонованого,
було отримано два головнi результати. По-перше, виведено явнi формули для елементiв обернених
матриць Мура–Пенроуза для двудiагональних матриць. По-друге, використовуючи цi формули,
побудовано скiнченний рекурсивний алгоритм, оптимальної обчислювальної складностi. Таким
чином, запропоновано варiант обчислення оберненої матрицi Мура–Пенроуза для двудiагональ-
них матриць, що не використовує сингулярний розклад.
Ключовi слова: псевдообернення Мура–Пенроуза, двудiагональна матриця, формула обер-
нення, рекурсивний алгоритм.

Матерiал надiйшов 28.06.2019