数理解析研究所講究録920 巻 1995 年 101-114

101

Abstract. A primal-dual framework of combinatorial relaxation algorithms is proposed for computing
the$\mathrm{h}$ighest degree of a minor of order$k$of a rational function matrix. The algorithm can be used forcomputing the index of nilpotency of a matrix pencil (or the index of the associated DAE). lt is a linear

algebraic version of the Hungarian method for the assignment problem. The proposed framework stands

in contrast with the previous combinatorial relaxation algorithm based on weighted matchings, and may
also be regarded as an extension of the Wolovich algorithm for$\mathrm{r}\mathrm{o}\mathrm{w}/\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{n}$properness. Several algorithms

are evaluated through computer experiments.
Key words.
$\mathrm{S}\mathrm{m}\mathrm{i}\mathrm{t}\mathrm{h}-\mathrm{M}\mathrm{C}\mathrm{M}\mathrm{i}|\mathrm{l}\mathrm{a}\mathrm{n}$combinatorial relaxation, degree of subdeterminant, index of DAE, Kronecker form,
form at infinity

11. 1 Introduction
Let $A(x)=(A_{\iota j}\langle x))$ be an $m\cross n$ rational function matrix with $A_{ij}(x)$ being a rational function in$x$with coefficients from the real number field R. This paper deals with algorithms for computing thehighest degree of a minor of order $k( \leq\min(m, n))$ of $A(x)$ :102$J]||I|=|J|=k$ },{deg det$AI,J[I,$(1)$\delta_{k}(A)=\max$where $A[I, J]$ denotes the submatrix of$A$with row-set

and column-set

$I$

$J$

and the degree of arational function $f(x)=\rho(x)/q(x)$ is defined by $\deg f(x)=\deg p(X)-\deg q(X)$ . By convention we put$\deg f(x)=-\infty$ if $f(x)=0$. This is one of the fundamental problems in combinatorial matrix theory [3]and has various applications in practice.
By computing$\delta_{k}\langle A$) $(k= 1,2, \cdots)$ we can obtain the

known also as the structure at infinity in the literature of control theory ($\mathrm{t}\mathrm{u}\mathrm{s}[12]$, Verghese-Kailath[22] . When $A(x)$ is a$)$$(k=1,2, \cdots)$

tions (DAEs)of order$\mathrm{r}\mathrm{e}\mathrm{g}\mathrm{u}\mathrm{l}\mathrm{a}\Gamma..\mathrm{p}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{i}.1$

determines the structural indices of its Kronecker form (seehence the index$\nu$form at infinity that is$\mathrm{S}\mathrm{m}\mathrm{i}\mathrm{t}\mathrm{h}- \mathrm{M}\mathrm{c}\mathrm{M}\mathrm{i}\mathrm{I}|\mathrm{a}\mathrm{n}$$n$$\mathrm{C}\mathrm{o}\mathrm{m}\mathrm{m}\mathrm{a}\mathrm{u}\mathrm{l}\mathrm{t}-\mathrm{D}\mathrm{i}\mathrm{o}\mathrm{n}[5]$, on the other hand,, Hau$\delta_{k}(A)$,Section 2.2]), and$\mathrm{M}\mathrm{u}\mathrm{r}\mathrm{o}\mathrm{t}\mathrm{a}[15$as well as the dynamical degree of reedom of the associated differential-algebraic equa$f$$(\mathrm{B}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{a}\mathrm{n}-\mathrm{c}_{\mathrm{a}}\mathrm{m}\mathrm{p}\mathrm{b}\mathrm{e}||-^{\mathrm{p}1}\mathrm{e}\mathrm{t}_{\mathrm{Z}}\mathrm{o}\mathrm{d}[2], \mathrm{G}\mathrm{e}\mathrm{a}\mathrm{r}[9,10], \mathrm{H}\mathrm{a}\mathrm{i}\Gamma \mathrm{e}\mathrm{r}^{-}\mathrm{W}\mathrm{a}\mathrm{n}\mathrm{n}\mathrm{e}\mathrm{r}[11])$. To be specific, the index is..

given by $\nu=\delta_{n-1}(A)-\delta_{\mathfrak{n}}(A)+1$ , whereas the dynamical degree of freedom is equal to$\delta_{n}(A)$. Thus the

problem of computing $\delta_{k}(A)(k=1,2, \cdots)$ has a fundamental engineering significance.
The present paper aims to propose a primal-d ual$\mathrm{f}_{\Gamma \mathrm{a}\mathrm{m}\mathrm{e}\mathrm{W}\mathrm{o}\mathrm{r}}\mathrm{k}$for computing the maximum degree of subdeterminants. Like theof combinatorial relaxation algorithms$\mathrm{H}$ungarian method for the assignmen$\mathrm{t}$problem $[1,20]$ , the algorithms keep a set of admissible potentials (feasible dual variables), which deter-

mines the estimate (dual objective value), and check its tightness (dual optimality) by solving a simpler(restricted primal) problem. If it turns out to be nontight, we obtain information useful for modifying thepotentials as well as the matrix. Thus the proposed .ramework can be interpreted as a linear algebraic$.\mathrm{f}$version of the Hungarian method. It stands in contrast with the previous combinatorial relaxation algorithm based on weighted-matchings, and may also be regarded as an extension of the Wolovich algorithm
for$\mathrm{r}\mathrm{o}\mathrm{w}/\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{n}$properness. Several algorithms are evaluated through computer experiments.

11.2 Combinatorial relaxation algorithms$\mathrm{l}\mathrm{n}$$\delta_{k}$this section, we present a new framework of$\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{b}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{t}_{0}\mathrm{r}\mathrm{i}\mathrm{a}1.$.

relaxation type algorithms for computing. This shares the approach of “combinatorial relaxation” with the previous$\mathrm{a}\downarrow \mathrm{g}_{\mathrm{o}\mathrm{r}}\mathrm{i}\mathrm{t}\mathrm{h}\dot{\mathrm{m}}$of$\mathrm{M}\mathrm{u}\mathrm{r}\mathrm{o}\mathrm{t}- \mathrm{a}[15]$,while featuring potentials (dual variables) rather than matchings in constructing the combinatorial estimate. The proposed framework may be thought of as a linear algebraic version of the primal-dual approachto the assignment problem that is well known in the context of combinatorial optimization $[1,20]$ .

103

11.2.1 Outline of the pr,o posed framework
For a rational function matrix $A(x)-$ , let $R$ and$C$denote the row-set and the column-set, respectively.

By potentials we mean integer numbers associated with rows or columns. With appropriately chosen
potentials, say, $\rho=(\rho_{i}|i\in R)$ and $\gamma=(\gamma_{j}|j\in C)$ , we compute an estimate

such that$\delta_{k}(A)\leq\wedge\delta_{k}(A;\rho, \gamma)$$\wedge\delta_{k}(A;\rho, \gamma)$

of$\delta_{k}(A)$. We then test for tightness (equality) in this inequality without computing. If the estimate turns out to be nontight, we modify the matrix $A(x)$ as well as the potentials.
$\mathrm{t}$$\delta_{k}(A)$The outline of the combinatorial relaxation algorithm to compute$\delta_{k}$for a fixed$k$is summarized asfollows:
Phase

$0$

: Choose appropriate$(\rho, \gamma)$.

Phase 1: Compute the combinatorial upper estimate
Phase 2: Test whether
Phase 3: Modify$\wedge\delta_{k}(A;\rho, \gamma)$

7$:=\gamma’$$A$$\wedge\delta_{k}(A\cdot\rho|’\gamma)=\delta_{k}(A)$$C_{1i}^{\cdot}$.

or not. If the equality holds, output

to another matrix $A’$ and$(\rho, \gamma)$to, namely $A’$ has the same true value assuch that $\delta_{k}(A’)=\delta_{k}(A)$ and$(\rho’, \gamma^{j})$$A$$\wedge\delta_{k}(A;\rho, \gamma)$

and stop.

$\wedge\delta_{k}(A’ ; \rho, \gamma^{l})’<$

and an improved estimate. Put $A:=A’,$ $\rho:=\rho’$ ,and go to Phase 1.

11.2.2 Phase 1–
The$\wedge\delta_{k}(A;\rho, \gamma)$$\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{b}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{i}\dot{\mathrm{a}}\downarrow$$\mathrm{C}\mathrm{o}\mathrm{m}..$estimate, denoted asbinatO-rial estimation$\wedge\delta_{k}(A;\rho, \gamma)$

in Section 11.2.1, is specified in this section. Define

by
’

:. .

”

$c_{ij}=\deg A_{ij}(X)$(1),

where by convention we put $c_{ij}=-\infty$ if $A_{ij}(x)=0$ . We$\mathrm{c}\mathrm{h}_{\mathrm{o}\mathrm{O}\dot{\mathrm{S}}}\mathrm{e}$a set of$\dot{\mathrm{p}}_{\mathrm{o}\mathrm{t}}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{a}$] $\mathrm{s}\rho 1=(\rho. |i\in R)$

and $\gamma=(\gamma_{j}|j\in C1$, which satisfies the condition$\rho_{\mathfrak{i}}+\gamma j\geq c_{i_{\mathrm{J}}}$

for $i\in R$ and $j\in C.$

Such potentials are sa..id to be admissible. Then for$I\subseteq R$

and$J\subseteq C$(2)with $|I|=|J|$ it holds that(3)$\deg\det A[I, J]\leq\rho(I)+\gamma(J)$ ,$\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}|\rho(I)=\sum_{i\in I}^{}\rho$

.

,.

and 7

$(J)= \sum_{J\in J}\gamma_{j}’$.

By considering the maximum possible value of the right-hand side, i.e., by introducing

104$\wedge\delta_{k}(A;\rho, 1)=\max\rho(I)+\max\gamma(J)$ ,$|J|=k$$|J|=k$we obviously have(4)$\rho(I)+\gamma(J)\leq\delta_{k}(A;\rho, \gamma\wedge)$for$I\subseteq R$

and$J\subseteq C$with $|I|=|J|=k$ . The inequalities (3) and (4) imply$\delta_{k}(A)\leq\delta_{k}(A;\rho, \gamma\wedge)$

which shows that$\wedge\delta_{k}(A;\rho, \gamma)$,(5)serves as an upper bound on the true value. Remember thatmust$(\rho, \gamma)$

be admissible in the sense of (2) before (5) can be true.

The estimate$\wedge\delta_{k}(A;\rho, \gamma)$

admits an alternative expression suitable for its efficient computation. Con-sider the following permutations for the row-set$\tau$$\sigma$: $\{1, \ldots, m\}arrow R$

: 1,
$\{$

$\ldots$, $n\}arrow C$

and the column-set $C$ of

$R$

$A$

:

such that $\rho\tau(1)\geq\rho\tau(2)\geq\cdots\geq\rho_{\tau(m)}$ ;

(6)

.

(7)

such that$\gamma_{\sigma(1)}\geq\gamma_{\sigma(2)}\geq\cdots\geq\gamma_{\sigma(n)}$Then it is easy to see that$\wedge\delta_{k}(A;\rho, \gamma)=\sum\rho_{\tau\langle:)}k+\sum\gamma\sigma(kj\rangle\cdot$$j=1$$.\cdot=1$

11.2.3 Phase 2–Test for tightness
This section provides a way to test whether$\delta_{k}(A)=\wedge\delta_{k}(A;\rho, \gamma)$

or not without computing$\delta_{k}(A)$directly.

Define$I^{t},$$J^{\mathrm{t}},\tilde{I},\overline{j}$as follows:$I^{t}=\{i\in R|\rho_{1}>\rho_{\mathcal{T}(}k+1)\}$,$\tilde{I}=\{i\in R|\rho i\geq\rho_{\mathcal{T}(k)}\}$$J^{\#}=\{j\in C|\gamma_{j}>\gamma_{\sigma(k+1})\}$,$\overline{J}=\{j\in C|\gamma_{j}\geq\gamma_{\sigma(k\rangle}\}$

;

(8)

.

(9)

Exceptionally, in case $k=m$ , we put $I\#=\overline{I}=R$ and in case $k=n$ , we put $J^{\#}=\tilde{J}=C$ . Then

attains the maximum subject to $|I|=k$ if and only if
to $|J|=k$ if and only if$J^{\mathrm{t}}\subseteq J\subseteq\tilde{J}$$J\#\subseteq J\subseteq\tilde{J}$. Define a constant matrix$I^{t}\subseteq I\subseteq\tilde{I}$,$\gamma(J)$$A^{t}$

byif$i\in\overline{I}$and$xarrow\infty$$A^{t}|.j=$otherwise.

$0$$\mathrm{Y}$$\#\not\in$

is maximized subject. Therefore (4) holds with equality if and only if$\lim x^{-\rho-\gamma}ijA_{j}.(x)$$A^{\#}=\langle A_{ij}^{\#})$. Similarly

11.2.1. (tightness) The following three conditions$(\mathrm{a})-(\mathrm{C})$are equivalent.

$\rho(I)$$j\in\overline{J}$;$I\#\subseteq I\subseteq\overline{I}$

and

105

(a)

$\delta_{k}(A)=\text{\’{o}}_{k}(A;\rho,\gamma)\wedge$(b) There exist$I\supseteq I^{\mathfrak{p}}$.
and$J\supseteq J\#$(c) The following four conditions

such that rank $A^{\rho}[I, J]=|I|=|J|=k$ .

$(\mathrm{r}\mathrm{l})-(\mathrm{r}4)$

are satisfied:(r1) ran $\mathrm{k}A\#[R, c]\geq k$ ;(r2) ran $\mathrm{k}A\#[I^{t}, c]=|I\#|_{1}$(r3) ran $\mathrm{k}A^{\#}[R, Jt]=|J\#|$ ;(r4) rank $A^{\#}[I^{\mathrm{t}}, J\#]\geq|I^{\mathfrak{g}}|+|J^{t}|-k$.

According to Proposition 11.2.1. above, the test for$\wedge\delta_{k}=\delta_{k}$in Phase 2 can be reduced to computing

the ranks of the four constant matrices. Actually, for (r1), (r2) and (r4), we use the Gaussian elimination
by elementary row transformations, and for (r3), we use the Gaussian elimination by elementary column

transformations.

11.2.4 Phase 3–Modification
When the estimate$\wedge\delta_{k}(A;\rho, \gamma)$

has turned out to be nontight, we modify the matrix as well as thepotentials. Recall that the maximum degree of subdeterminants of a fixed order is invariant under biproper
equivalence transformations. That is,$A’(x)=U(x)A(x)V(x)$(10)

satisfies $\delta_{k}(A’)=\delta_{k}(A)$ if both $U(x)$ and $V(x)$ are biproper. We use this type of transformation in
the modification of the matrix $A(x)$ .

Let

$U$

be a nonsingutar constant matrix of order $m$ such that
$U_{h}.\cdot=0$

if$\rho_{h}>\rho:$,and define a transformation matrix $U(x)$ by$U(x)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(x;\rho)\cdot U$. diag $(x;-\rho)$ ,where diag $(x;\rho)$ designates the diagonal matrix diagmatrix. Similarly, using a nonsingular constant matrix$V_{jh}=0$

if

$V$

(11)

$(x^{\rho_{1}} , x^{\rho_{2}}, \cdots , x^{\rho_{m}})$. Then $U(x)$ is a biproperthat satisfies the condition:$\gamma_{h}>\gamma j$,

we define a transformation matrix $V(x)$ by$V(x)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\langle_{X;-}\gamma$

which is also a biproper matrix.

)$\cdot$V. diag$(x;\gamma)$

,

(12)

106

We now explain how to construct such a pair of constant matrices

which of the four conditions$(\mathrm{r}\mathrm{l})-(\mathrm{r}4)$

is violated by the constant matrix

elimination with column pivoting, applied to the matrix

that SA$\mathfrak{p}$$\rho|$

we obtain a nonsingular matrix’

has ( $m$ –rank$|H|=|\tilde{I}|-k+1$

and$A^{t}$$V$

efficiently according to(see Proposition 11.2.1.

$(\mathrm{c})$

$A\#$

$S$

with rows rearranged by$A\#$of order $m$ such that$\tau$

according to$S_{\tau \mathrm{t}^{1)_{T}\langle i)}}=0$

if$\mathit{1}<|$and (SA ) $[H, C]=O$ . Construct the constant matrix

where $I_{|R-}H|$ denotes the unit matrix. Let$U$

as follows:

pivoting, applied to the matrix$A^{t}[I^{\mathrm{t}}, C]$

we obtain a nonsingular matrix–rank

row vector in$A^{t_{[I^{t},c}}]$$\mathrm{l}\mathrm{n}$

$n$

has (

$n$

.

$|J^{\mathrm{t}}|$$|I\#|$

such that$U$$\tau$

according to the magnitude of$s_{\tau(\iota)}\tau(:)=0$

if $l<i$ and that $S\cdot A\#[I\#, C]$

as in (13). Let$V$

be the unit matrix of order$n$.

$A\#[R, J^{\mathrm{t}}]$–rank $A^{t}[R,$$J^{1}]$

with columns rearranged by$S$of ordersuch that$|J^{\mathrm{t}}|$$\sigma$

according to the magnitude$S_{\sigma(\mathrm{j})\sigma}(t)=0$

if 1 $<j$ and that) zero column vectors. Construct the constant matrix
$U$

of

(14)

be the unit matrix of order $m$ .

When (r4) is violated:$\mathrm{l}\mathrm{n}$the course of checking the condition (r4) by means of the Gaussian row

elimination with column pivoting, applied to the matrix
to the magnitude ofand that$V$in the case of (r2), namely$V[C, H]=S[C, H],$ $V[H, C-H]=O,$ $V[C-H, C-H]=I_{|C-H|}$ .

$U$,

be the singleton set corresponding to a zero

$H\subseteq I\#$in a similar (transposed) way as we constructedLet$\rho_{i}$checking the condition (r3) by the Gaussian column elimination with row, we obtain a nonsingular matrix$A\#[R, J^{t}]\cdot S$

order

of order. Then constructpivoting, applied to the matrix$\gamma j$$S$with rows rearranged by) zero row vectors. Let$S\cdot A^{1}[I^{\mathrm{t}}, C]$

When (r3) is violated:ofbe the unit matrix of order

$V$

(13)

In checking the condition (r2) by the Gaussian row elimination with column

When (r2) is violated:$|I\#|$

and) zero vectors among its row vectors. Let $H$ be a subset of \‘I such that

$U[H, R]=S[H, R],$ $U[R-H, H]=O,$ $U[R-H, R-H]=I_{|R-}H|$ ,

has (

).

In the course of checking the condition (r1) by means of the Gaussian row

When (r1) is violated:the magnitude of$U$$\rho:$$S\cdot A^{t}[I^{t}, J\#]$, we obtain a nonsingular matrix

has ($|I\#|$-rank $A\#[I\#,$

to $(k-|J^{t}|+1)$ zero row vectors in

way as in (13). Let$V$$J\#]$$A^{\mathrm{t}}[I^{1}, J^{t}]$$S$

with rows rearranged byof order$|I\#|$

such that$I\#$. Then construct the constant matrix

be the unit matrix of order$n$according$S_{\mathcal{T}(\downarrow)}\tau(i)=0$) zero row vectors. Let $H$ be a subset of$S\cdot A^{t}[I\#, J^{\mathrm{t}}]$$\tau$

if$\mathit{1}<i$corresponding$U$in the same.

After modifying the matrix $A(x)$ to $A’(x)$ by (10) with $U(x)$ and $V(x)$ thus constructed, the set of
potentials is updated as follows. Putthat$c_{\mathrm{j}}’.\cdot=\deg A_{i}’\mathrm{j}(x)$similarly to (1). It follows from (10), (11) and (12)107diag $(x;-\rho)\cdot A^{J}(x)$ . diag $(x;-\gamma)=U$ . diag $(x;-\rho)\cdot A(X)$ . diag $(x;-\gamma)\cdot V$,
which, together with (2), implies$\rho$.

$+\gamma_{j}\geq c_{ij}’$

for $i\in R$

and $j\in C.$(15)

We update the set of potentials to improve the estimate exploiting the fact that (15) above holds withstrict inequality for some $i\in R$ and $j\in C$ . Roughly speaking, we reduce the potentials of such rows orcolumns without violating the admissibility. More specific description is given below according to which

condition of the four is violated.
When (r1) is violated:$\max\{\max(j\in^{\overline{r}}C_{j}’.\cdot-\gamma_{g}\cdot),j\in(c-\overline{r})\max C\prime ij-\gamma_{\sigma(k)}\}$$\rho_{i}’$

$=$

$=$

$i\in\tilde{I}$;(16)otherwise;$\rho_{i}$$\gamma_{j}’$

if

$\{$

$\max_{1\in R}(C_{j}^{J}.\cdot-\rho_{1}’.)$

.

(17)

When (r2) is violated:$\max\{\max(cj\in\tilde{J}’|\mathrm{J} -\gamma_{\mathrm{J}}),\max c^{J}J\in(c-\overline,)|\mathrm{J} -\gamma_{\sigma(k)}\}$$\rho_{i}’$$=$I$\iota_{;}$

otherwise;$\rho_{i}$$=$$i\in$(18).

$\gamma_{j}’$

ir

$\{$

$\max(c_{i}’i\in Rj-\rho_{i}’)$

.

(19)

When (r3) is violated:$\max\{.\cdot\max(c’\in f\dot{.}j-\rho_{i}),\max c_{i}^{J}i\in(R-\overline{I})j-\rho\overline{.}(k)\}$$\gamma_{j}’$

$=$

(20)

otherwise;$\gamma j$$\rho_{\mathrm{i}}’$$=$if $j\in J\#$ ;$\{$$\max(cj\in CJi\mathrm{j}-\gamma_{j}’)$

.

(21)

When (r4) is violated:. $-\beta$$\rho|$

if $i\in H$ ;$\gamma_{j}$

if $j\in J\#$ ;(22)$\gamma_{j}’=\{$$\rho_{i}’=\{$$\rho_{i}$otherwise,$\gamma_{j}+\beta$

otherwise,

where$\beta=\min\{\rho_{\mathcal{T}}(|I\#|)-\rho_{\tau(1|}I\#+1)$ ’$\gamma_{\sigma(|J|)}\#-\gamma_{\sigma(\mathrm{I}J|)}\#+1$,$\min_{h\in H.j\in J\#}(\rho h+\gamma_{j}-c_{hj}’)\}$

(23)

.

Exceptionally, the first term is excluded from the minimization when $I^{t}=R$ and the second is

excluded when $J^{t}=C$ .
As to the admissibility of the updated set of potentials, the following lemma holds.

108

11.2.2. The updated set of potentials$\ovalbox{\tt\small REJECT}\#$

That is,$\rho_{i}’+\gamma_{j}’\geq \mathrm{C}_{1j}’.$$(\rho’, \gamma)’$

is admissible with respect to the modified matrix $A’$ .

for $i\in R$ and $j\in C$ .

The following theorem shows that the estimate is improved by the modification descri bed above.
11.2.3. The following inequality holds:$\mathrm{E}\mathrm{E}$$\gamma)\wedge$$\wedge\delta_{k}(A^{lJ\prime} ; \rho, \gamma)<\delta_{k(}A;\rho,$

.

(24)

11.2.5 Algorithm description
We have so far fixed the size

obtaining the]]$\mathrm{S}\mathrm{m}\mathrm{i}\mathrm{t}\mathrm{h}-\mathrm{M}\mathrm{c}\mathrm{M}\mathrm{i}$$\mathrm{a}\mathrm{n}$$k$of a submatrix. It is necessary to computefor $k=1,2,$$\delta_{k}$$\ldots$

,

$r$

for

form at infinity.

Before describing the whole algorithm, we should remark how we can detect$\delta_{k}=-\infty$. For a rational

fu nction$f(x)= \frac{a_{1}x^{d_{1}}+O_{2^{X^{d_{2}d}}}+\cdots+a_{l^{X}}p}{q(x)}$$(d_{1}>d_{2}>\cdots>dp)$ ,we define min-deg $f(x)\equiv d_{l}-\deg q(x)$ with respect to its expression. Put
while $d_{\max}=., \max_{1j}\deg Aij(x)$ . Ir$\delta_{k}=-\infty$

when$\delta_{k}$$\wedge\delta_{k}(A;\rho, \gamma)<k\cdot d_{\mathrm{m}\iota \mathrm{n}}$is finite, it must satisfy
for some admissible

We now summarize the algorithm for computing
not necessary to input$r$$\delta_{k}$$\delta_{k}\geq k\cdot d_{\mathrm{m}1\mathrm{n}}$$(\rho, \gamma)$$d_{\mathrm{m}\mathrm{i}\mathfrak{n}}= \min_{!^{j}}$,min-deg $A_{j}.(x)$. Hence we can conclude that.

for $k=1,2,$, since it is computed in the procedure.

$\ldots$, ( $=$ rank $A$ )$r$

as follows. It is

Algorithm for computing$\delta_{k}(A)$$(k=1,2, \ldots)$

Step

$0$

:

Compute$d_{\min}$. Set $k:=1$ . Find a set of initial potentials$(\rho, \gamma)$

which satisfies$\rho,$$+\gamma_{j}\geq c_{1j}$for $i\in R$ and $j\in C$ .

Step 1:

Compute$\wedge\delta_{k}:=\sum\rho_{\tau(}:)k+\sum\gamma_{\sigma}\langle kj)$

$i=1$

where

.. .

$\geq\gamma_{\sigma(\mathfrak{n})}$Step 2:$\tau$. If

and$\sigma$$j=1$

are permutations such that$\wedge\delta_{k}<k\cdot d_{\min}$,$\rho_{\tau(1)}\geq\rho_{f}(2)\geq\cdots\geq\rho_{\tau\langle m)}$, then halt (at this point $(k-1)$ equals the rank

$r$

and

).

Put
;

$:=\{i\in R|\rho_{i}>\rho\tau(k+1)\}$,$\tilde{I}:=\{i\in R|\rho:\geq\rho_{\tau(k)}\}$$J^{t}:=\{j\in C|\gamma_{j}>\gamma_{\sigma \mathrm{t}^{k+1})}\}$,$\tilde{J}:=\{j\in C|\gamma_{j}\geq\gamma\sigma(k)\}$$I^{\#}$;

where exceptionally [ $\#=R$ and $J^{t}=C$ in case $k=m$ and $k=n$ , respectively. Put$\gamma_{\sigma(1)}\geq\gamma_{\sigma(2)}\geq$109$\lim x^{-\rho:-\gamma}jA|.j(_{X)}$$A^{\#}.\cdot j:=\{$

If the four conditions

if$i\in\overline{I}$and$j\in\tilde{J}$;$xarrow\infty$$0$

otherwise.

$(\mathrm{r}\mathrm{l})-(\mathrm{r}4)$

in Proposition 11.2.1. are satisfied, then output$\wedge\delta_{k}$

.

$1f$

at leas$\mathrm{t}$one condition of the four is violated, go to Step 3. If $k= \min(m, n)$ , then halt. Otherwise $k:=k+1$
and go to Step 1.

Step 3:Modify the matrix $A(x)$ and the potentials as described in (16) $-(22)$ , according to which of$(\mathrm{r}1)-(\Gamma 4)$

is violated. Go to Step 1.

In the algorithm described above, we may start with an arbitrary set of admissible potentials. It is
desirable, however, to adopt reasonably small integers as the initial potentials.

One of the natural choices is to put$\rho_{i}=\max c_{j}j\in C|$

for $i\in R$ ,$\gamma_{j}=0$

for $j\in C.$

In this case, $\gamma_{j}=0$ for $j\in C$ throughout the computation. Hencewith the row degrees. Since$J^{t}=\emptyset$(25)$\rho_{i}$for $i\in R$ always coincide

and $\tilde{J}=C$ , the four rank conditions in the test for tightness (seeProposition 11.2.1. ) are reduced to a couple of conditions (r1) and (r2). This algorithm will be called

the row-degree method, and is an extension of the well-known algorithm for$\delta_{m}(A)$of an $m\cross n$ row-fullrank matrix $A(x)$ based on the concept of row properness (Wolovich [23]).

Another natural choice is, of course, the symmetric (transposed) version of the above, i.e.,$\rho_{i}=0$

for $i\in R$ ,.

$\gamma_{j}=.\max_{1\in R}c_{j}$for $j\in C$ .

(26)

This algorithm, however, does not behave in a symmetric manner to the row-degree method. Thisis because the symmetry between rows and columns is broken in the way of modification. In fact, any
condition of the four may possibly be violated. Thus this algorithm is not the “column-degree method,”
and the primal-dual framework is not a simple extension of the Wolovich algorithm.

11.3 Experimental evaluation

11.3.1 Method
We make experiments on Sun SPARC Station 10 to evaluate the three combinatorial relaxation type

110

algorithms, the primal-dual method$(\mathrm{R}\mathrm{D})$$(\mathrm{P}\mathrm{D})$with initial potentials given by (26), the row-degree method

and the matching method (M). All of these are implemented inlanguage. For comparison we also$\mathrm{C}$realize the elimination method (E) in Maple (a computer algebraic system).

Our program restricts input matrices to Laurant polynomial matrices of which each entry is a polynomial in

and $1/x$ .

$x$We provide two types of sample problems, regular pencils (P) and band matrices (B). A regular pencil$xE-F$ used here is generated randomly in the following way:$E_{:j}=\{$

with probability

9$\frac{1}{2}-\not\in$for $g=\pm 1,$$\pm 2,$

for $h=\pm 1,$$\pm 2,$$\ldots$$,$$\ldots$$,$$\pm K$;$\pm h’$;with probability ,$0$$\rho$

with probability $*1-2$$h$$F_{j}.\cdot=\{$with probability

$0$

$q$

.

We put $h\cdot=1,$ $\rho=0.0625$ and $q=0.5$ .

On the other hand, we provide band matrices in the following way. First we put$X=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(x;d)$

where

$d$

.

,

$d=(d_{1}, d_{2}, \ldots d_{\mathfrak{n}})$,is generated randomly as follows:$d_{:}=\ell$

with probability$\frac{1}{2L+1}$

for $\ell=0,$$\pm 1,$ $\pm 2,$$\ldots$$\pm L$$,$.

With a lower bidiagonal matrix$\frac{1}{x}1$

$P=($

$.1.$

.

$)$

..

$\cdot$$\frac{1}{x}$1and an upper bidiagonal matrix

1$. \frac{1}{x}..\cdot$

.

$Q=($

$)$

$.1^{\cdot}$$\frac{1}{x}1$

we construct a pentadiagonal (band) matrix $A(x)=P^{2}XQ^{2}$ . We put $L=10$ . Obviously, both
$Q$

is

$P$

and

are biproper, and hence the Smith-McMillan form at infinity of $A(x)$ is the same as that of $X$ , which$\mathrm{e}\mathrm{a}s$ily obtained by sorting$(d_{1}, \ldots , d_{\mathfrak{n}})$. Therefore, we know$\delta_{k}(A)$in advance and we use this type of

111

matrices as sample data to encertain that our program returns correct answers.
We measure$\mathrm{C}$PU time and the number of modifications (for$\mathrm{R}\mathrm{D}$$\mathrm{P}\mathrm{D},$, and M) to evaluate eachalgorithm.

11.3.2 Results
Before we report experimental results, we should note that the time complexity is $\vee(rmn\max(m, n))$if $(d_{\max}-d_{\min})$ is bounded by a constant. In fact, we have $(d_{\max}-d_{\min})\leq 1$ with the pencils and$(d_{\max}-d_{\min})\leq 2L+2$$\not\equiv 1$

with the band matrices. Therefore the time complexity is V $(rmn \max(m, n))$ .

Statistics over 10 pencils of each size.
$\cap \mathrm{f}$$\mathrm{m}\cap r\mathrm{i};\mathrm{f}\mathrm{l}p\mathrm{a}*\mathrm{i}\mathrm{n}\mathrm{n}\epsilon$$\mathrm{e}\mathrm{h}\cap \mathrm{w}$$*\mathrm{h}$a$\mathrm{v}\mathrm{n}1\cdot$$\mathfrak{m}\mathrm{a}\mathrm{n}$$\mathrm{n}$CPU time in seconds and the number, aa$\cap \mathrm{f}$1$\cap$s-$\mathrm{i}\mathrm{t}_{\theta \mathrm{P}\mathrm{f}\mathrm{i}}\mathrm{t}.\mathrm{i}\cap \mathfrak{n}$

[ $V$

$=$

1

$\mathfrak{n}$$=$$\cap\cap 69^{\not\in}.\mathrm{i}$

$n$

$=$

Tables 1 and 2 show the results for the pencils and the band matrices, respectively. Algorithm$(\mathrm{P}\mathrm{D})$large

failed for $n=128$ because of a memory shortage. The experiments with algorithm (E) for$n$were given up because it would have taken many days.

It is observed in our experiment that for$n\cross n$matrix$A$, approximately V $(n^{3.4})$ time is neededby each of combinatorial relaxation algorithm, while the worst case time complexity is$\vee(n^{4})$sincerank $A=n$ both in problem type (P) and in (B). On the other hand, the elimination method

takes approximately
functions is$\vee(n^{s_{)}}$.

$\vee(n^{4})$time in our experiment, whereas the number of operations on rational

112$\Leftrightarrow’$$.\mathrm{q}_{l}-\mathrm{f}\mathrm{i}\prime \mathrm{i}\mathrm{s}’-\mathrm{i}r\mathrm{C}\cap\vee 0r1\cap \mathrm{h}$a$\mathrm{n}’ \mathrm{i}\mathrm{m}\mathrm{n}*r\mathrm{i}$,os$\cap \mathrm{f}\mathrm{o}$fi$r\mathrm{h}\propto \mathrm{i}\sigma \mathrm{o}$$\cap \mathrm{D}$I$\mathrm{T}*\cdot \mathrm{m}\mathrm{o}$;$\mathrm{n}\mathrm{e}\mathrm{o}r\cap \mathrm{n}$$A$ $\epsilon\circ \mathrm{n}’ \mathrm{I}*\mathrm{h}$a$\mathrm{n}\cdot,$$\mathrm{m}\mathrm{k}L\mathrm{r}-\mathrm{f}$In the above experiment, no distinction in computing time is recognized among algorithms$(\mathrm{P}\mathrm{D}),$$(\mathrm{R}\mathrm{D})$, and (M). Since algorithm$(\mathrm{R}\mathrm{D})$lacks symmetry, it is natural that this algorithm

takes longer time than the other algorithms for a matrix with a special structure. For example,$A(x)=$
let us consider

(1)

$F3$

CPU time in seconds and the number of modiflcations for the matrix $A(x)$ of (1) with
$n=16$.

The results for the matrix with $n=16$ are shown on Table 3.

113

11.4 Conclusion
In this paper, we have proposed a primal-dual framework of combinatorial relaxation algorithms for computing the maximum degree of subdeterminants. Furthermore, we have realized
combinatorial relaxation algorithms in computer program and made experiments. It turns outthat combinatorial relaxation algorithms are effective for computing the maximum degree of sub-determinants and the index of nilpotency.

参考文献$[1]\mathrm{R}$. K. AHLJA, T. L. MAGNANTI, J. B. ORLIN, Network Flows –Theory, Algorithms andApplications, Prentice Hall, 1993.
$[2]\mathrm{K}$. E. BRENAN, S. L. CAMPBELL, A$l\mathfrak{i}\mathrm{D}$L. R. PETZOLD, Numerical Solution of Initial-Value

Problems in Differential-Algebraic Equations, North-Holland, 1989.
[3] R. A.

$\mathrm{B}\mathrm{R}\mathrm{L}\mathrm{A}\mathrm{L}\mathrm{D}\mathrm{l}\mathrm{A}\ltimes \mathrm{D}$H. J. RYSER, Combinatorial Matrix Theory, Cambridge University Press,

1991.
[4] P. BL JAKIEW[C’Z, Maximum Weighted Matching for High Index Differential Algebraic Equa-tions, Doctor Thesis, Delft University of Technology, 1994.
$[5]\mathrm{C}$. COMMAU$\mathrm{L}\mathrm{T}\mathrm{A}\uparrow \mathfrak{i}\mathrm{D}$J.-M. DION, Structure at infinity of linear multivariable systems: a geo-

metric approach, IEEE Trans. Automat. Control, AC-27 (1982), pp. 693-696.
$[6]\mathrm{C}$. COMMAU$\mathrm{L}\mathrm{T}$, J.-M. DION A$b\mathfrak{i}\mathrm{D}$A. PEREZ, Disturbance rejection for structured systems,

IEEE Trans. Automat. Control, AC-36 (1991), pp. 884-887.
$[7]1$.

$\mathrm{D}\mathrm{t}\mathrm{F}\mathrm{F}\mathrm{A}\aleph^{\vee}\mathrm{D}$C. W. GEAR, Computing the structural index, SIAM J. Algebraic Discrete

Methods, 7 (1986), pp. 594-603.
$[8]\mathrm{R}$. GANI$\mathrm{A}1\backslash \mathrm{D}$I. T. CAMERON, Modelling for dynamic simulation of chemical processes: theindex problem, Chem. Engin. Sci., 47 (1992), pp. 1311-1315.
$[9]\mathrm{C}$. W. GEAR, Differential-algebraic equation index transformations, SIAM J. Sci. Stat. Com-put., 9 (1988), pp. 39-47.
[$101\mathrm{c}$. W. GEAR, Differential algebraic equations, indices, and integral algebraic equations, SIAM

J. Numer. Anal., 27 (1990), pp. 1527-1534.
$[11]\mathrm{E}$

. HA

1991.

$1\mathrm{R}\mathrm{E}\mathrm{R}\mathrm{A}\ltimes \mathrm{D}\mathrm{G}^{1}$.

$\backslash \iota^{\gamma}\mathrm{A}\ltimes|\backslash \mathrm{E}\mathrm{R}$, Solving Ordinary Differential Equations II, Springer-Verlag,

114. L. J. HALTLS, The forrnal Laplace transform for smooth linear systems, in Mathematical. G. Marchesini and S. K. Mitter), Lecture Notes in Economics andSystems Theory ($[12]\mathrm{M}$$\mathrm{e}\mathrm{d}\mathrm{s}$Mathematical Systems, 131, Springer-Verlag, 1976, pp. 29-47.
$[13]\mathrm{K}$. MuROTA, Systems Analysis by Graphs and Matroids –Stru ctural Solvability and Control-

lability, Algorithms and Combinatorics 3, Springer-Verlag, 1987.

.

$[14]\mathrm{K}$$\mathrm{M}\mathrm{t}$ROTA, Computing Puiseux-series solutions to deterrninantal equations via combinato rial

relaxation, SIAM J. Comput., 19 (1990), pp. 1132-1161.

.

$[15]\mathrm{K}$

ROTA, Combinatorial relaxation algorithm$\mathrm{M}\mathrm{L}$for the maximum degree of subdeterminants:

Computing Smith-McMillan forrn at infinity and structural indices in Kronecker form, Appl.

Algebra Engin. Comm. Comput., to appear.
[16] K. ML ROTA, Computing the degree of determinants via combinatorial relaxation, SIAM J.

Comput., 24 (1995), No. 4, to appear.
MLROTA AND J. W.

$[17]\mathrm{K}.$$\mathrm{V}\mathrm{A}\aleph$DER$\dagger \mathrm{V}\mathrm{o}\iota^{\vee}\mathrm{o}\mathrm{E}$, Structure at infinity of structured descriptor systems

and its applications, SIAM J. Control Optim., 29 (1991), pp. 878-894.
$[18]\mathrm{M}$. NEWMAK, Integral Mat rices, Academic Press, 1972.

. C.

$[19]\mathrm{C}$$\mathrm{P}\mathrm{A}\aleph \mathrm{T}\mathrm{E}\mathrm{L}1\mathrm{D}\mathrm{E}\mathrm{S}$, The consistent initializationof differential-algebraic systems, SIAM

J. Sci. Stat. Comput., 9 (1988), pp. 213-231.
. H.

$[20]\mathrm{C}$$\mathrm{p}_{\mathrm{A}\mathrm{p}\mathrm{A}\mathrm{D}\mathrm{l}\mathrm{M}}\mathrm{I}\mathrm{T}\mathrm{R}\mathrm{I}\mathrm{O}\{_{-}\mathrm{A}\ltimes \mathrm{D}$K. STEIGLITZ, Combinato rial Optimization: Algorithms and Com-

plexity, Prentice Hall, 1982.
$[21]\mathrm{J}$.

$\mathrm{U}\ltimes \mathrm{G}\mathrm{A}\mathrm{R}$, A. KR\"ONER AND W. MARQLARDT, Structural analysis of differential-algebraic

equation systems: Theory and application, Comput. Chem. Engin., to appear.
$[22]\mathrm{G}$. C. VERGHESE AND T. KAILATH, Rational matrix structure, IEEE Trans. Automat. Control,

AC-26 (1981), pp. 434-439.
$[23]\mathrm{W}$$[24]\mathrm{J}$

. A.

. W.

$\mathrm{W}\mathrm{o}\mathrm{L}\mathrm{O}\backslash j$

1CH, Linear Multivariable Systems, Springer-Verlag, 1974.

$\mathrm{V}\mathrm{A}\ltimes \mathrm{D}\mathrm{E}\mathrm{R}$WOUDE, On the structure at infinity of a structured system, Linear AlgebraAppl., 148 (1991), pp. 145-169.