Inversion of Analytic Matrix Functions That are Singular at the Origin.
In this paper we study the inversion of an analytic matrix valued function A(z). This problem can also be viewed as an analytic perturbation of the matrix A0=A(0). We are mainly interested in the case where A0 is singular but A(z) has an inverse in some punctured disc around z=0. It is known that A-1(z) can be expanded as a Laurent series at the origin. The main purpose of this paper is to provide efficient computational procedures for the coefficients of this series. We demonstrate that the proposed algorithms are computationally superior to symbolic algebra when the order of the pole is small.
Introduction
Let fA k g k=0;1;::: ' R n\Thetan be a sequence of matrices that de-nes the analytic matrix valued
function
The above series is assumed to converge in some non-empty neighbourhood of z = 0. We will
also say that A(z) is an analytic perturbation of the matrix A Assume the inverse
matrices A \Gamma1 (z) exist in some (possibly punctured) disc centred at In particular, we
are primarily interested in the case where A 0 is singular. In this case it is known that A \Gamma1 (z)
can be expanded as a Laurent series in the form
A
and s is a natural number, known as the order of the pole at z = 0. The main
purpose of this paper is to provide eOEcient computational procedures for the Laurent series
coeOEcients As one can see from the following literature review, few computational
methods have been considered in the past.
This work was supported in part by Australian Research Council Grant #A49532206.
y INRIA Sophia Antipolis, 2004 route des Lucioles, B.P.93, 06902, Sophia Antipolis Cedex, France, e-mail:
k.avrachenkov@sophia.inria.fr
z Department of Statistics, The Hebrew University, 91905 Jerusalem, Israel and Department of Economet-
rics, The University of Sydney, Sydney, NSW 2006, Australia, e-mail: haviv@mscc.huji.ac.il
x CIAM, School of Mathematics, The University of South Australia, The Levels, SA 5095, Australia, e-mail:
The inversion of nearly singular operator valued functions was probably -rst studied in
the paper by Keldysh [22]. In that paper he studied the case of a polynomial perturbation
are compact operators on Hilbert space. In particular, he showed that
the principal part of the Laurent series expansion for the inverse operator A \Gamma1 (z) can be
given in terms of generalized Jordan chains. The generalized Jordan chains were initially
developed in the context of matrix and operator polynomials (see [13, 26, 30] and numerous
references therein). However, the concept can be easily generalized to the case of an analytic
perturbation (1).
Following Gohberg and Sigal [15] and Gohberg and Rodman [14], we say that the vectors
Jordan chain of the perturbed matrix A(z) at
for each 0 - k - r \Gamma 1. Note that ' 0 is an eigenvector of the unperturbed matrix A 0
corresponding to the zero eigenvalue. The number r is called the length of the Jordan chain
and ' 0 is the initial vector. Let f' (j)
j=1 be a system of linearly independent eigenvectors,
which span the null space of A 0 . Then one can construct Jordan chains initializing at each
of the eigenvectors ' (j)
0 . This generalized Jordan set plays a crucial role in the analysis of
analytic matrix valued functions A(z).
Gantmacher [11] analysed the polynomial matrix (3) by using the canonical Smith form.
Vishik and Lyusternik [37] studied the case of a linear perturbation
showed that one can express A \Gamma1 (z) as a Laurent series as long as A(z) is invertible in some
punctured neighbourhood of the origin. In addition, an undetermined coeOEcient method for
the calculation of Laurent series terms was given in [37]. Langenhop [25] showed that the
coeOEcients of the regular part of the Laurent series for the inverse of a linear perturbation
form a geometric sequence. The proof of this fact was re-ned later in Schweitzer [33, 34] and
Schweitzer and Stewart [35]. In particular, the paper [35] proposed a method for computing
the Laurent series coeOEcients. However the method of [35] cannot be applied (at least imme-
diately) to the general case of an analytic perturbation. Many authors have obtained existence
results for operator valued analytic and meromorphic functions [3, 15, 23, 27, 29, 36]. In par-
ticular, Gohberg and Sigal [15], used a local Smith form to elaborate on the structure of the
principal part of the Laurent series in terms of generalized Jordan chains. Recently, Gohberg,
Kaashoek and Van Schagen [12] have re-ned the results of [15]. Furthermore, Bart, Kaashoek
and Lay [5] used their results on the stability of the null and range spaces [4] to prove the
existence of meromorphic relative inverses of -nite meromorphic operator valued functions.
The ordinary inverse operator is a particular case of the relative inverse. For the applications
of the inversion of analytic matrix functions see for example [8, 9, 20, 23, 24, 28, 31, 32, 36].
Howlett [20] provided a computational procedure for the Laurent series coeOEcients based
on a sequence of row and column operations on the coeOEcients of the original power series
(1). Howlett used the rank test of Sain and Massey [32] to determine s, the order of the pole.
He also showed that the coeOEcients of the Laurent series satisfy a -nite linear recurrence
relation in the case of a polynomial perturbation. The method of [20] can be considered as a
starting point for our research. The algebraic reduction technique which is used in the present
paper was introduced by Haviv and Ritov [17, 18] in the special case of stochastic matrices.
Haviv, Ritov and Rothblum [19] also applied this approach to the perturbation analysis of
semi-simple eigenvalues.
In this paper we provide three related methods for computing the coeOEcients of the Laurent
series (2). The -rst method uses generalized inverse matrices to solve a set of linear
equations and extends the work in [17] and [20]. The other two methods use results that
appear in [2, 17, 18, 19] and are based on a reduction technique [6, 10, 21, 23]. All three
methods depend in a fundamental way on equating coeOEcients for various powers of z. By
substituting the series (1) and (2) into the identity A(z)A I and collecting coeOE-
cients of the same power of z, one obtains the following system which we will refer to as the
fundamental equations:
A similar system can written when considering the identity A I but of course
the set of fundamental equations (4:0); suOEcient. Finally, for matrix operators
each in-nite system of linear equations uniquely determines the coeOEcients of the Laurent
series (2). This fact has been noted in [3, 20, 23, 37, 36].
Main results
Let us de-ne the following augmented matrix A (t) 2 R (t+1)n\Theta(t+1)n
A (t) =6 6 6 6 6 4
A t A
and prove the following basic lemma.
s be the order of the pole at the origin for the inverse function A \Gamma1 (z). Any
eigenvector \Phi 2 R (s+1)n of A (s) corresponding to the zero eigenvalue possesses the property
that its -rst n elements are zero.
Proof: Suppose on the contrary that there exists an eigenvector \Phi 2 R (s+1)n such that
A
and not all of its -rst n entries are zero. Then, partition the vector \Phi into s blocks and
rewrite (5) in the form
with ' 0 6= 0. This means that we have found a generalized Jordan chain of length s + 1.
However, from the results of Gohberg and Sigal [15], we conclude that the maximal length
of a generalized Jordan chain of A(z) at z = 0 is s. Hence, we came to a contradiction and,
consequently,
direct proof of Lemma 1 is given in Appendix 1.
vectors \Phi 2 R (s+j+1)n in the null space of the augmented matrix A (s+j) , j - 0,
possess the property that the -rst (j + 1)n elements are zero.
The following theorem provides a theoretical basis for the recursive solution of the in-nite
system of fundamental equations (4).
Theorem 1 Each coeOEcient X k , k - 0 is uniquely determined by the previous coeOEcients
and the set of s fundamental equations
Proof: It is obvious that the sequence of Laurent series coeOEcients fX i
i=0 is a solution
to the fundamental equations (4). Suppose the coeOEcients X i ,
determined. Next we show that the set of fundamental equations (4.k)-(4.k+s) uniquely
determines the next coeOEcient X k . Indeed, suppose there exists another solution ~
are both solutions, we can write
A (s)6 4
~
~
and
A (s)6 4
where the matrix J i is de-ned as follows:
ae
and where ~
X k+s are any particular solutions of the nonhomogenous linear system
(4.k)-(4.k+s). Note that (6) and (7) have identical righthand sides. Of course, the dioeerence
between these two righthand sides, [ ~
is in the right null space of
A Invoking Lemma 1, the -rst n rows of [ ~
are hence zero. In
other words, ~
which proves the theorem.Using the above theoretical background, in the next section we provide three recursive
computational schemes which are based on the generalized inverses and on a reduction tech-
nique. The reduction technique is based on the following result. A weaker version of this
result was utilized in [17] and in [19].
Theorem 2 Let fC k g t
suppose that the
system of t equations
is feasible. Then the general solution is given by
(R
where C y
0 is the Moore-Penrose generalized inverse of C 0 and Q 2 R m\Thetap is any matrix whose
columns form a basis for the right null space of C 0 . Furthermore, the sequence of matrices
solves a reduced -nite set of t matrix equations
where the matrices D k 2 R p\Thetap and S k 2 R p\Thetan , are computed by the following
recursion. Set U
Then,
where M 2 R p\Thetam is any matrix whose rows form a basis for the left null space of C 0 .
Proof: The general solution to the matrix equation (7.0) can be written in the form
arbitrary matrix.
In order for the equation
to be feasible, we need that the right hand side R belongs to R(C 0
is
where the rows of M form a basis for N(C T
Substituting expression (13) for the general
solution into the above feasibility condition, one -nds that W 0 satis-es the equation
which can be rewritten as
Thus we have obtained the -rst reduced fundamental equation (9.0) with
Next we observe that the general solution of equation (7.1) is represented by
the formula
(R
with . Moving on and applying the feasibility condition to equation (7.2), we
obtain
and again the substitution of expressions (13) and (14) into the above condition yields
(R
which is rearranged to give
The last equation is the reduced equation (9.1) with
. Note that this equation imposes restrictions on W 1 as well as on
By proceeding in the same way, we eventually obtain the complete system of equations
with coeOEcients given by formulas (11) and (12) each of which can be proved by induction
in a straightforward way.Remark 3 In the above theorem it is important to observe that the reduced system has the
same form as the original but the number of matrix equations is decreased by one and the
coeOEcients are reduced in size to matrices in R p\Thetap , where p is the dimension of N(C 0 ) or,
equivalently, the number of redundant equations de-ned by the coeOEcient C 0 .
In the next section we use this reduction process to solve the system of fundamental
equations. Note that the reduction process can be employed to solve any appropriate -nite
subset of the fundamental equations.
3 Solution methods
In this section we discuss three methods for solving the fundamental equations. The -rst
method is based on the direct application of Moore-Penrose generalized inverses. The second
method involves the replacement of the original system of the fundamental equations by
a system of equations with a reduced dimension. In the third method we show that the
reduction process can be applied recursively to reduce the problem to a non-singular system.
Since all methods depend to some extent on the prior knowledge of s, we begin by discussing
a procedure for the determination of s. A special procedure for determining this order for the
case where the matrices A(z) are stochastic and the perturbation series is -nite is given in [16].
It is based on combinatorial properties (actually, network representation) of the processes and
hence it is a stable procedure. However, as will be seen in Section 3.4, it is possible to use
the third method without prior knowledge of s. Actually, the full reduction version of our
procedure determines s as well. Of course, as in any computational method which is used to
determine indices which have discrete values, using our procedures in order to compute the
order of singularity might lack stability.
3.1 The determination of the order of the pole
The rank test on the matrix A (t) proposed by Sain and Massey in [32] is likely to be the
most eoeective procedure for determining the value of s. The calculation of rank is essentially
equivalent to the reduction of A (t) to a row echelon normal form and it can be argued that
row operations can be used successively in order to calculate the rank of A (0) ,A (1) ,A (2)
and -nd the minimum value of t for which rankA (t\Gamma1) +n. This minimum value of t equals s,
the order of the pole. Note that previous row operations for reducing A (t\Gamma1) to row echelon
form are replicated in the reduction of A (t) and do not need to be repeated. For example,
if a certain combination of row operations reduces A 0 to row echelon form, then the same
operations are used again as part of the reduction of
to row echelon form.
3.2 Basic generalized inverse method
In this section we obtain a recursive formula for the Laurent series coeOEcients X k , k - 0 by
using the Moore-Penrose generalized inverse of the augmented matrix A (s) .
y be the Moore-Penrose generalized inverse of A (s) and de-ne the matrices
G
G
G
0s
G
Furthemore, we would like to note that in fact we use only the -rst n rows of the generalized
namely, [G
Proposition 1 The coeOEcients of the Laurent series (2) can be calculated by the following
recursive formula
s
G
0s and the matrix J i is de-ned by
ae
Proof: According to Theorem 1, once the coeOEcients X i , are determined, the
next coeOEcient X k can be obtained from the (4.k)-(4.k+s) fundamental equations.
A (s)6 4
The general solution to the above system is given in the form6 6 6 4
~
~
G
0s
G
1s
G
where the -rst block of matrix \Phi is equal to zero according to Lemma 1. Thus, we immediately
obtain the recursive expression (15). In particular, applying the same arguments as above to
the -rst s we obtain that
0s .Note that the matrices J j+k in the expression (15) disappear when the regular coeOEcients
are computed.
Remark 4 The formula (15) is a generalization of the recursive formula for the case where
A 0 is invertible. In this case,
while initializing with
Remark 5 Probably from the computational point of view it is better not to compute the
generalized inverse G (s) beforehand, but rather to -nd the SVD or LU decomposition of A
and then use these decompositions for solving the fundamental equations (3:k)-(3:k + s). This
is the standard approach for solving linear systems with various righthand sides.
3.3 The one step reduction process
In this section we describe an alternative scheme that can be used in the case where it is
relatively easy to compute the bases for the right and for the left null spaces of A 0 . Speci-cally,
be the dimension of the null space of A 0 , let Q 2 IR n\Thetap be a matrix whose
columns form a basis for the right null space of A 0 and let M 2 IR p\Thetan be a matrix whose p
rows form a basis for the left null space of A 0 . Of course, although
possible, we are interested in the singular case where p - 1.
Again, as before, we suppose that the coeOEcients X i ,
determined. Then, by Theorem 1, the next coeOEcient X k is the unique solution to the
subsystem of fundamental equations
The above system is like the one given in (9) with C and with R
Therefore, we can apply the reduction process described
in Theorem 2. This results in the system
where the coeOEcients D i and S i , can be calculated by the recursive formulae
(11) and (12).
Remark 6 Note that in many practical applications p is much less than n and hence the
above system (17) with D i 2 IR p\Thetap is much smaller than the original system (16).
Now we have two options. We can either apply the reduction technique again (see the
next subsection for more details) or we can solve the reduced system directly by using the
generalized inverse approach. In the latter case, we de-ne
and
0t
Then, by carrying out a similar computation to the one presented in the proof of Proposition 1,
we obtain
Once W 0 is determined it is possible to obtain X k from the formula
Furthermore, substituting for S i , 0 - i - s\Gamma1, from (12) and changing the order of summation
gives
A y
s
Note that by convention the sum disappears when the lower limit is greater than the upper
limit. Now, substituting R
into the expression (18), we
obtain the explicit recursive formula for the Laurent series coeOEcients
A y
s
A (J k+j \Gamma
for all k - 1. In particular, the coeOEcient of the -rst singular term in (2) can be given by the
3.4 The complete reduction process
As was pointed out in the previous section, the reduced system has essentially the same
structure as the original one and hence one can apply again the reduction step described in
Theorem 2. Note that each time the reduction step is carried out, the number of matrix
equations is reduced by one. Therefore one can perform up to s reduction steps. We now
outline how these steps can be executed. We start by introducing the sequence of reduced
systems. The fundamental matrix equations for the l-th reduction step are
A (l)
A (l)
A (l)
s\Gammal X (l)
s\Gammal
one gets the original system of fundamental equations and with gets the
reduced system for the -rst reduction step described in the previous subsection. Initializing
with R (0)
I and with A (0)
s, the matrices A (l)
and R (l)
for each reduction step 1 - l - s, can be computed successively by a
recursion similar to (11) and (12). In general we have
U (l)
A (l\Gamma1)
A (l)
R (l)
U (l)
where Q (l) and M (l) are the basis matrices for the right and left null spaces respectively of
the matrix A (l\Gamma1)
0 and where A (l\Gamma1)y
0 is the Moore-Penrose generalized inverse of A (l\Gamma1)
. After
s reduction steps, one gets the -nal system of reduced equations
A
is a unique solution to the subsystem of fundamental equations (4.0)-(4.s) and
Theorem 2 states the equivalence of the l-th and (l 1)-st systems of reduced equations, the
system (22) possesses a unique solution, and hence matrix A
0 is invertible. Thus,
The original solution X
0 can be now retrieved by the backwards recursive relationship
Now by taking R (0)
one gets the algorithm for computing
the Laurent series coeOEcients recursive formulae similar to (15) and (19)
can be obtained, but they are quite complicated in the general case.
The order s of the pole can also be obtained from the reduction process by continuing the
process until A (l)
becomes non-singular. The number of reduction steps equals the order of
the pole. Note also that the sequence of matrices A (l)
can be computed irrespectively
of the right hand sides. Once s is determined, one can compute R (l)
Computational complexity and comparison with symbolic algebr

In this section we compare the computational complexity of the one-step-reduction process
when applied to compute X 0 with the complexity of symbolic algebra. In particular, we show
that the former comes with a reduced complexity in the case where the pole has a relatively
small order. The computational complexity of the other two procedures can be determined
similarly.
To compute the coeOEcients D i , of the reduced fundamental system (17),
one needs to perform O(s 2 n 3 ) operations. The total number of reduced equations is sp
(recall that p is the dimension of the null space of A 0 ). Hence, the computational complexity
for determining X 0 by the one-step-reduction process is O(maxfs 2 g). The Laurent
series (2) in general, and the coeOEcient X 0 in particular, can also be computed by using
symbolic algebra. This for example can be executed by MATLAB symbolic toolbox and is
done as follows. Since X 0 is uniquely determined by the -rst s equations
one needs to do in order to compute X 0 is it to invert symbolically the
following matrix polynomial
Symbolic computations here mean performing operations, such as multiplication and division,
over the -eld of rational functions (and not over the -eld of the reals). In particular, if the
degrees of numerators and of denominators of rational functions do not exceed q, then each
operation (multiplication or division) which is performed in the -eld of rational functions
translates into qlog(q) operations in the -eld of real numbers [1]. Note that during the
symbolic inversion of the polynomial matrix (25), the degree of rational functions does not
exceed sn. The latter fact follows from Cramer's rule. Thus, the complexity of the symbolic
inversion of (25) equals O(n 3 ) \Theta log(sn)). As a result, one gets a matrix
A \Gamma1 (z) whose elements are rational functions of z. The elements of the matrix X 0 can then be
immediately calculated by dividing the leading coe-cients of the numerator and denominator.
Finally, one can see that if s !! n and p !! n, which is typically the case, then our method
comes with a reduced computational burden.
Concluding remarks
In this paper we have shown that the Laurent series for the inversion of an analytic matrix
valued function can be computed by solving a system of fundamental linear equations.
Furthermore, we demonstrated that the system of fundamental equations can be solved recur-
sively. In particular, the coeOEcient X k is determined by the previous coeOEcients X
and the next s is the order of the pole. We suggest three
basic methods, one without any reduction (see (15)), one with a single reduction step (see
and (20)), and one using a complete reduction process with s steps (see (23) and (24)).
Of course, an intermediate process with the number of reductions between 1 and s could be
used too. We note that when the complete reduction process is used the order of the pole
can be determined through the execution of the algorithm. When s !! n and p !! n, the
proposed algorithms by far outperform the method based on symbolic algebra.

Acknowledgement

The authors are grateful to Prof. Jerzy A. Filar for his helpful advice. Also the authors would
like to thank anonymous referees for their valuable suggestions and for directing us to some
existing literature.
Apendix 1: Another proof of Lemma 1
A direct proof of Lemma 1 can be carried out using augmented matrices. Speci-cally, de-ne
are the coeOEcients of the Laurent series (2). Then it follows from
the fundamental systems (4) and (5) that the augmented matrices A (t) and X (t) satisfy the
relationship
where the augmented matrix E (t) 2 R (t+1)n\Theta(t+1)n is de-ned by setting
p;q=0 where
n\Thetan and
ae
I for
s:
Now, as before, let \Phi 2 R (s+1)n satisfy the equation
A
If we multiply equation (27) from the left by X reduces to
The vector E (s) \Phi has ' 0 as the (s 1)-st block, which gives the required result.
Apendix 2: A Numerical example
Let us consider the matrix valued function
where 2. Construct the augmented matrices
and note that which is the dimension of the original
coeOEcients A 0 and A 1 . Therefore, according to the test of Sain and Massey [32], the Laurent
expansion for A \Gamma1 (z) has a simple pole. Alternatively, we can compute a basis for
which in this particular example consists of only one vector
\Theta
The -rst three zero elements in q (1) con-rm that the Laurent series has a simple pole. Next
we compute the generalized inverse of A (1) given by
G (1)
1=3 \Gamma5=12 \Gamma1=12 1=8 1=8 \Gamma1=8
1=3 \Gamma5=12 \Gamma1=12 \Gamma3=8 \Gamma3=8 3=8
Consequently,
\Gamma3 \Gamma3
Alternatively, we know that X 0 is uniquely determined by the fundamental equations
After one reduction step these equations reduce to
where
\Theta

and
\Theta
Hence,
\Theta

and

\Gamma3 \Gamma3
The latter expression is identical with (28) and coincides with the one computed by expanding
A \Gamma1 (z) with the help of the MATLAB symbolic toolbox. Note that even for this three
dimensional example the direct symbolic calculation of the Laurent series takes a relatively
long time.



--R

The design and analysis of computer algorithms
The fundamental matrix of singularly perturbed Markov chains
Meromorphic operator valued functions.
Stability properties of
Relative inverses of meromorphic operator functions and associated holomorphic projection functions

Generalized Inverses of Linear Transformation
iSingular systems of dioeerential equationsj
iSingular systems of dioeerential equations IIj
A reduction process for perturbed Markov chains
The theory of matrices
On the local theory of regular analytic matrix functions
Matrix Polynomials
Analytic matrix functions with prescribed local data
An operator generalization of the logarithmic residue theorem and the theorem of Rouch
Mean passage times and nearly uncoupled Markov chains
Series expansions for
Matrix Anal.
iTaylor expansions of eigenvalues of perturbed matrices with applications to spectal radii of nonnegative matrices
Input retrieval in
Perturbation theory for linear operators
On the characteristic values and characteristic functions of certain classes of non-selfadjoint equations
Mathematical foundations of the state lumping of large systems
iInversion of lambda-matrices and application to the theory of linear vibrationsj
The Laurent expansion for a nearly singular matrix
Introduction to the spectral theory of polynomial operator pencils
iSpectral properties of a polynomial op- eratorj
Theory of Suboptimal Decisions

An introduction to operator polynomials
The Laurent expansion of a generalized resolvent with some applications
Invertibility of linear time invariant dynamical systems
The Laurent expansion for a nearly singular pencil
Perturbation series expansions for nearly completely-decomposable Markov chains
The Laurent expansion of pencils that are singular at the origin
Theory of branching of solutions of non-linear equations
The solution of some perturbation problems in the case of matrices and self-adjoint and non-self-adjoint dioeerential equations
--TR

--CTR
Jerzy A. Filar, Controlled Markov chains, graphs, and Hamiltonicity, Foundations and Trends in Stochastic Systems, v.1 n.2, p.77-162, January 2006
matrix inversion;matrix valued functions;analytic perturbation;laurent series