A Fast Algorithm of Moore-Penrose Inverse for the Loewner-Type Matrix

Abstract

In this paper, we present a fast algorithm of Moore-Penrose inverse for 𝑚×𝑛 Loewner-type matrix with full column rank by forming a special block matrix and studying its inverse. Its computation complexity is 𝑂(𝑚𝑛)+𝑂(𝑛²), but it is 𝑂(𝑚𝑛²)+𝑂(𝑛³) by using 𝐋⁺=(𝐋^T𝐋)⁻¹𝐋^T.

1. Introduction

Loewner matrix was first studied by Loewner in 1934 in [1]. He studied the relations of various Loewner matrices via the characteristic of the monotone matrix function and the problem of rational interpolation at that time. Since then, more further studies had been carried out by many scientists in [2–5], such as Belevitch, Donoghue Jr., Fiedler, Chen and Zhang. From their mentioned works, we can find various properties of Loewner matrix and its application in the rational interpolation. In 1984, Vavrín presented a fast algorithm for the inverse of Loewner matrix in [6], then the fast algorithm for the Loewner system was got accordingly. Rost and Vavrín presented a fast algorithm for the system whose coefficient is a Loewner-Vandermonde matrix from 1995 to 1996 in [7, 8]. Lu gave a fast triangular factorization algorithm for the symmetrical Loewner-type matrix in 2003 in [9]. Xu et al. and so forth gave a fast triangular factorization algorithm for the inverse of symmetrical Loewner-type matrix in 2003 in [10]. In 2009, Tong et al. gave a fast algorithm of the Moore-Penrose inverse for 𝑚×𝑛 symmetrical Loewner-type matrix [11]. In this paper, we generalize the fast algorithm of symmetrical Loewner-type matrix to the Loewner-type matrix. The theory and computation of generalized inverse matrix arise in various applications such as optimization theory, control theory, computation mathematics, and mathematical statistics. Therefore, there are very important theoretical and practical significance when we study the fast algorithm of Moore-Penrose inverse for the Loewner-type matrix.

This paper is organized as follows: in Section 1, we present some preliminaries results. The fast algorithm of Moore-Penrose inverse for Loewner-type matrix is driven in Section 2. We give some numerical examples to illustrate the fast algorithm obtained in Section 3.

2. Preliminaries

A Loewner matrix 𝐋₁ is a matrix of the form 𝐋₁=(𝑙_𝑖𝑗)_𝑚×𝑛=((𝜉_𝑖−𝜂_𝑗)/(𝛼_𝑖−𝛽_𝑗))^𝑚,𝑛_𝑖,𝑗=1, where 𝛼_𝑖,𝛽_𝑗,𝜉_𝑖,and𝜂_𝑗 (𝑖=1,2,…,𝑚;𝑗=1,2,…,𝑛) are given numbers and 𝛼_𝑖≠𝛽_𝑗. A Loewner-type matrix 𝐋 is a matrix which satisfies 𝐃𝐋−𝐋𝐃=^𝑙_𝑗=1𝐩^(𝑗)𝐪^(𝑗)T, (2.1) where 𝐃=diag(𝛼₁,𝛼₂,…,𝛼_𝑚), 𝐃=diag(𝛽₁,𝛽₂,…,𝛽_𝑛), 𝐩^(𝑗)=(𝑝₁^(𝑗),𝑝₂^(𝑗),…,𝑝_𝑚^(𝑗))^T, and 𝐪^(𝑗)=(𝑞₁^(𝑗),𝑞₂^(𝑗),…,𝑞_𝑛^(𝑗))^T. Loewner matrix 𝐋₁ is a special case of Loewner-type matrix, and it satisfies 𝐃𝐋₁−𝐋₁𝐃=(𝜉₁,…,𝜉_𝑚)^T(1,…,1)−(1,…,1)^T(𝜂₁,…,𝜂_𝑛).

Let 𝑛 be the rank of 𝑚×𝑛 Loewner-type matrix 𝐋. Forming an 𝑚+𝑛 matrix 𝐌=−𝐈_𝑚𝐋𝐋^T𝟎, (2.2) we obtain𝐌⁻¹=𝐋𝐋^𝑇𝐋⁻¹𝐋^T−𝐈_𝑚𝐋𝐋^T𝐋⁻¹𝐋^𝑇𝐋⁻¹𝐋^T𝐋^T𝐋⁻¹=𝐋𝐋⁺−𝐈_𝑚𝐋^T+𝐋⁺𝐋^T𝐋⁻¹. (2.3)

It is obvious that we may obtain the Moore-Penrose inverse 𝐋⁺ of Loewner-type matrix 𝐋 by (2.3) if we can get 𝐌⁻¹.

Now let us begin to find 𝐌⁻¹.

The following result will be useful for getting 𝐌⁻¹.

Lemma 2.1 (see [12]). Let all the leading principal submatrices 𝐒_𝑖(𝑖=1,2,…,𝑛) of matrix 𝐒=(𝑠_𝑖𝑗)_𝑛×𝑛 whose order is 𝑛 be invertible. Linear system 𝐒𝐱=𝐛 be given, where 𝐛=(𝑏₁,𝑏₂,…,𝑏_𝑛)^T. Note that 𝐛_𝑖=(𝑏₁,𝑏₂,…,𝑏_𝑖)^T, and let 𝐱_𝑖=(𝑥_𝑖1,𝑥_𝑖2,…,𝑥_𝑖𝑖)^𝑇, 𝐮_𝑖=(𝑢_1𝑖,𝑢_2𝑖,…,𝑢_𝑖𝑖)^T be solution vectors of linear systems 𝐒_𝑖𝐱_𝑖=𝐛_𝑖, 𝐒_𝑖𝐮_𝑖=𝐞_𝑖^(𝑖) differently, then 𝐱_𝑖=𝐱_𝑖−10+𝜎_𝑖𝐮_𝑖, (2.4) where 𝜎_𝑖=𝑏_𝑖−∑^𝑖−1_𝑗=1𝑠_𝑖𝑗𝑥_{𝑖−1,𝑗}.

3. Fast Algorithm of the Moore-Penrose Inverse for Loewner-Type Matrix

By using (2.1), we know that 𝐌 satisfies 𝐃𝐃𝐃𝐃𝐌−𝐌=^𝑙_𝑗=1𝐩^(𝑗)𝟎_𝑛𝟎^T_𝑚𝐪^(𝑗)T−𝟎_𝑚𝐪^(𝑗)𝐩^(𝑗)T𝟎^T_𝑛. (3.1)

Let all the leading principal submatrices 𝐌_𝑖??(𝑖=1,2,…,𝑚+𝑛) of 𝐌 be invertible. If 𝑖>𝑚, by virtue of (3.1), we have𝐃𝐃_𝑖−𝑚𝐌_𝑖−𝐌_𝑖𝐃𝐃_𝑖−𝑚=^𝑙_𝑗=1𝐩^(𝑗)𝟎_𝑖−𝑚𝟎^T_𝑚𝐪^(𝑗)T_𝑖−𝑚−𝟎_𝑚𝐪^(𝑗)_𝑖−𝑚𝐩^(𝑗)T𝟎^T_𝑖−𝑚. (3.2)

Let 𝐠_𝑖^(𝑗)=(𝑔^(𝑗)_𝑖1,𝑔^(𝑗)_𝑖2,…,𝑔^(𝑗)_𝑖𝑖)^T, 𝐡_𝑖^(𝑗)=(ℎ^(𝑗)_𝑖1,ℎ^(𝑗)_𝑖2,…,ℎ^(𝑗)_𝑖𝑖)^T be solution vectors of linear systems 𝐌_𝑖𝐠_𝑖^(𝑗)=𝐩^(𝑗)𝟎_𝑖−𝑚,𝐌_𝑖𝐡_𝑖^(𝑗)=𝟎_𝑚𝐪^(𝑗)_𝑖−𝑚(𝑖=𝑚,𝑚+1,…,𝑚+𝑛) (3.3) differently, then using 𝐌_𝑚=−𝐈_𝑚, we obtain 𝐠_𝑚^(𝑗)=−𝐩^(𝑗), 𝐡_𝑚^(𝑗)=𝟎_𝑚. By virtue of Lemma 2.1 we have 𝐠_𝑖^(𝑗)=𝐠^(𝑗)_𝑖−10+𝜎_𝑖^(𝑗)𝐮_𝑖,𝐡_𝑖^(𝑗)=𝐡^(𝑗)_𝑖−10+𝜏_𝑖^(𝑗)𝐮_𝑖, (3.4) where 𝜎_𝑖^(𝑗)=−∑^𝑚_𝑘=1𝑙_{𝑘,𝑖−𝑚}𝑔^(𝑗)_{𝑖−1,𝑘}, 𝜏_𝑖^(𝑗)=𝑞^(𝑗)_𝑖−𝑚−∑^𝑚_𝑘=1𝑙_{𝑘,𝑖−𝑚}ℎ^(𝑗)_{𝑖−1,𝑘}, and 𝐮_𝑖=(𝑢_1𝑖,𝑢_2𝑖,…,𝑢_𝑖𝑖)^T is a solution vector of 𝐌_𝑖𝐮_𝑖=𝐞_𝑖^(𝑖). Multiplying (3.2) by 𝐌_𝑖⁻¹ on the left and on the right differently, we obtain𝐌_𝑖⁻¹𝐃𝐃_𝑖−𝑚−𝐃𝐃_𝑖−𝑚𝐌_𝑖⁻¹=^𝑙_𝑗=1𝐠_𝑖^(𝑗)𝐡^(𝑗)T_𝑖−𝐡_𝑖^(𝑗)𝐠^(𝑗)T_𝑖. (3.5) Multiplying (3.5) by 𝐞_𝑖^(𝑖) and noting that 𝐃_{𝐃𝑖−𝑚}𝐞_𝑖^(𝑖)=𝛽_𝑖−𝑚𝐞_𝑖^(𝑖), we have 𝛽_𝑖−𝑚𝐮_𝑖−𝐃𝐃_𝑖−𝑚𝐮_𝑖=^𝑙_𝑗=1ℎ^(𝑗)_𝑖𝑖𝐠_𝑖^(𝑗)−𝑔^(𝑗)_𝑖𝑖𝐡_𝑖^(𝑗). (3.6) Substituting (3.4) in (3.6) and noting that 𝑔^(𝑗)_𝑖𝑖=𝜎_𝑖^(𝑗)𝑢_𝑖𝑖, ℎ^(𝑗)_𝑖𝑖=𝜏_𝑖^(𝑗)𝑢_𝑖𝑖, we have𝛽_𝑖−𝑚𝐮_𝑖−𝐃𝐃_𝑖−𝑚𝐮_𝑖=𝑢^𝑙_{𝑖𝑖𝑗=1}𝜏_𝑖^(𝑗)𝐠^(𝑗)_𝑖−10−𝜎_𝑖^(𝑗)𝐡^(𝑗)_𝑖−10, (3.7) and hence, 𝑢_𝑘𝑖=𝑢_𝑖𝑖∑^𝑙_𝑗=1𝜏_𝑖^(𝑗)𝑔^(𝑗)_{𝑖−1,𝑘}−𝜎_𝑖^(𝑗)ℎ^(𝑗)_{𝑖−1,𝑘}𝛽_𝑖−𝑚−𝛼_𝑘(),𝑢𝑘=1,2,…,𝑚_𝑖𝑖∑^𝑙_𝑗=1𝜏_𝑖^(𝑗)𝑔^(𝑗)_{𝑖−1,𝑘}−𝜎_𝑖^(𝑗)ℎ^(𝑗)_{𝑖−1,𝑘}𝛽_𝑖−𝑚−𝛽_𝑘−𝑚(𝑘=𝑚+1,𝑚+2,…,𝑖−1). (3.8) Now, let us look for 𝑢_𝑖𝑖. Choosing the 𝑖th of 𝐌_𝑖𝐮_𝑖=𝐞_𝑖^(𝑖), namely ^𝑚_𝑘=1𝑙_{𝑘,𝑖−𝑚}𝑢_𝑘𝑖, (3.9) and using (3.8), we have 𝑢_𝑖𝑖=1∑^𝑚_𝑘=1𝑙_{𝑘,𝑖−𝑚}/𝛽_𝑖−𝑚−𝛼_𝑘∑^𝑙_𝑗=1𝜏_𝑖^(𝑗)𝑔^(𝑗)_{𝑖−1,𝑘}−𝜎_𝑖^(𝑗)ℎ^(𝑗)_{𝑖−1,𝑘}. (3.10)

Note that 𝐌⁻¹=(𝐯₁,𝐯₂,…,𝐯_𝑚+𝑛), where 𝐯_𝑖=(𝑣_1𝑖,𝑣_2𝑖,…,𝑣_{𝑚+𝑛,𝑖})^T is the 𝑖th column of 𝐌⁻¹, then 𝐯_𝑚+𝑛=𝐮_𝑚+𝑛. (3.11) Letting 𝑖=𝑚+𝑛 in (3.5) and multiplying it by 𝐞_𝑘^(𝑚+𝑛)(𝑘>𝑚) on the right, we have𝛽_𝑘−𝑚𝐯_𝑘−𝐃𝐃𝐯_𝑘=^𝑙_𝑗=1ℎ^(𝑗)_{𝑚+𝑛,𝑘}𝐠^(𝑗)_𝑚+𝑛−𝑔^(𝑗)_{𝑚+𝑛,𝑘}𝐡^(𝑗)_𝑚+𝑛. (3.12)

So, we obtain the fast algorithm of Moore-Penrose inverse 𝐋⁺ of Loewner-type matrix 𝐋 by (2.3), (3.4), and (3.8)~(3.12).

Algorithm
Step 1. 𝐠_𝑚^(𝑗)=−𝐩^(𝑗),𝐡_𝑚^(𝑗)=𝟎_𝑚𝑗=1,2,…,𝑙, (3.13)
For 𝑖=𝑚+1,𝑚+2,…,𝑚+𝑛, 𝜎_𝑖^(𝑗)=−^𝑚_𝑘=1𝑙_{𝑘,𝑖−𝑚}𝑔^(𝑗)_{𝑖−1,𝑘},𝜏_𝑖^(𝑗)=𝑞^(𝑗)_𝑖−𝑚−^𝑚_𝑘=1𝑙_{𝑘,𝑖−𝑚}ℎ^(𝑗)_{𝑖−1,𝑘},𝑡𝑗=1,2,…,𝑙_𝑖𝑘=^𝑙_𝑗=1𝜏_𝑖^(𝑗)𝑔^(𝑗)_{𝑖−1,𝑘}−𝜎_𝑖^(𝑗)ℎ^(𝑗)_{𝑖−1,𝑘}(),𝜆𝑘=1,2,…,𝑖−1_𝑖=^𝑚_𝑘=1𝑙_{𝑘,𝑖−𝑚}𝛽_𝑖−𝑚−𝛼_𝑘𝑡_𝑖𝑘,𝑢_𝑖𝑖=1𝜆_𝑖,𝑢_𝑘𝑖=𝑢_𝑖𝑖𝛽_𝑖−𝑚−𝛼_𝑘𝑡_𝑖𝑘(),𝑢𝑘=1,2,…,𝑚_𝑘𝑖=𝑢_𝑖𝑖𝛽_𝑖−𝑚−𝛽_𝑘−𝑚𝑡_𝑖𝑘(),𝐠𝑘=𝑚+1,𝑚+2,…,𝑖−1_𝑖(𝑗)=𝐠^(𝑗)_𝑖−10+𝜎_𝑖^(𝑗)𝐮_𝑖,𝐡_𝑖^(𝑗)=𝐡^(𝑗)_𝑖−10+𝜏_𝑖^(𝑗)𝐮_𝑖𝑗=1,2,…,𝑙. (3.14)Step 2. 𝑣_{𝑖,𝑚+𝑛}=𝑢_{𝑖,𝑚+𝑛}(𝑖=1,2,…,𝑚). (3.15)
For 𝑘=𝑚+1,…𝑚+𝑛−1, 𝑣_𝑖𝑘=1𝛽_𝑘−𝑚−𝛼_𝑖^𝑙_𝑗=1𝑔^(𝑗)_{𝑚+𝑛,𝑖}ℎ^(𝑗)_{𝑚+𝑛,𝑘}−ℎ^(𝑗)_{𝑚+𝑛,𝑖}𝑔^(𝑗)_{𝑚+𝑛,𝑘}(𝑖=1,…,𝑚), (3.16) then 𝐋⁺=𝑣_1,𝑚+1⋯𝑣_𝑚,𝑚+1𝑣⋮⋮_1,𝑚+𝑛⋯𝑣_{𝑚,𝑚+𝑛}. (3.17)
The algorithm requires (7𝑙−3)𝑚𝑛+𝑙𝑛²+𝑂(𝑚) multiplication and division operations and (5𝑙−3)𝑚𝑛+(𝑙/2)𝑛²+𝑂(𝑚) addition and subtraction operations, and the computation complexity is 𝑂(𝑚𝑛)+𝑂(𝑛²), but it is 𝑂(𝑚𝑛²)+𝑂(𝑛³) by using 𝐋⁺=(𝐋^T𝐋)⁻¹𝐋^T.

4. Numerical Examples

We get the Moore-Penrose inverse matrix of Loewner-type matrix with Fortran program in computer, and the results of the partial numerical examples are given as follows in Table 1 (the error is measured by the 2-norm of vector, and the time is measured by second).

Table 1 

Example 1. Consider 𝛼_𝑖=(𝑖−1)𝜋𝑚−𝑛+1,𝛽_𝑗=𝑗+1𝜋𝑚+𝑛−1,𝜉_𝑖=(−1)^𝑖(𝑖−𝑚𝑛),𝜂_𝑗=𝑗^𝑗−𝑚.𝑝⁽¹⁾(𝑖)=𝜉_𝑖2,𝑝⁽²⁾(𝑖)=1,𝑝⁽³⁾(𝑖)=𝜉_𝑖2,𝑝⁽⁴⁾(𝑖)𝑞=1;₍₁₎𝑗=1,𝑞⁽²⁾𝑗=−𝜂_𝑗,𝑞⁽³⁾𝑗=1,𝑞⁽⁴⁾𝑗=−𝜂_𝑗.𝑖=1,2,…,𝑚;𝑗=1,2,…,𝑛. (4.1)

From above examples and many more examples not given, we find that the stability of the fast algorithm is very good, and it needs shorter time than that of 𝐋⁺=(𝐋^𝑇𝐋)⁻¹𝐋^𝑇.

Acknowledgment

This work is being supported by the National Natural Science Foundation Grants of China (no. 60974082).