text:
\section{Introduction}
An $l$-th order \emph{differential-algebraic equation} (DAE) for $\funcdoms{x}{\setR}{\setR^n}$ is a differential equation in the form of
\begin{align} \label{def:dae}
  \app[\big]{F}{t, \app{x}{t}, \app{\dot{x}}{t}, \ldots, \app{x^{(l)}}{t}} = 0,
\end{align}
where $\funcdoms{F}{\setR \times \setR^n \times \cdots \times \setR^n}{\setR^n}$ is a sufficiently smooth function.
DAEs have aspects of both ordinary differential equations (ODEs) $\app{\dot{x}}{t} = \app{\phi}{t, \app{x}{t}}$ and algebraic equations $\app{G}{t, \app{x}{t}} = 0$.
DAEs are widely used for modeling of dynamical systems such as mechanical systems, electrical circuits and chemical reaction plants.

A DAE~\eqref{def:dae} can be rewritten as a first order DAE
\begin{align} \label{def:1st_dae}
  \app{F}{t, \app{x}{t}, \app{\dot{x}}{t}} = 0
\end{align}
by replacing higher order derivatives of $x$ with newly introduced variables.
The difficulty in numerically solving the DAE~\eqref{def:1st_dae} is measured by its \emph{differentiation index}~\cite{Brenan1996}, which is defined as the minimum nonnegative integer $\nu$ such that the system of equations
\begin{align*}
  \app{F}{t, \app{x}{t}, \app{\dot{x}}{t}} = 0, \quad
  \dif{}{t} \app{F}{t, \app{x}{t}, \app{\dot{x}}{t}} = 0, \quad
  \ldots, \quad
  \dif{^\nu}{t^\nu} \app{F}{t, \app{x}{t}, \app{\dot{x}}{t}} = 0
\end{align*}
can determine $\dot{x}$ as a continuous function of $t$ and $x$.
In other words, $\nu$ is the number of times one has to differentiate the DAE~\eqref{def:1st_dae} to get an ODE.
Intuitively, the differentiation index represents how far the DAE is from ODEs.

A common approach for solving a high $\prn{\geq 2}$ index DAE is to convert it into a low $\prn{\leq 1}$ index DAE.
This process is called an \emph{index reduction}, which is important for accurate simulation of dynamical systems.
Most of existing simulation software packages for dynamical systems, such as Dymola, OpenModelica, MapleSim and Simulink, are equipped with the index reduction algorithm given by Mattsson--S\"{o}derlind~\cite{Mattsson1993} (MS-algorithm).
The MS-algorithm uses Pantelides' method~\cite{Pantelides1988} as preprocessing.
Pantelides' method constructs a bipartite graph from structural information of a given DAE, and solves an assignment problem on the bipartite graph efficiently.
The MS-algorithm then differentiates equations in the DAE with the aid of the information obtained by Pantelides' method, and replaces some derivatives with dummy variables.
The MS-algorithm returns a sparse DAE if a given DAE is sparse, and thus the algorithm can be applied to large scale DAEs.

Pantelides' method, however, does not work even for the following simple DAE
\begin{align*}
  \begin{cases}
    \dot{x}_1 + \dot{x}_2 + x_3 = 0, \\
    \dot{x}_1 + \dot{x}_2 \phantom{{}+x_3} = 0, \\
    \phantom{x_1+{}} x_2 + \dot{x}_3 = 0. \\
  \end{cases}
\end{align*}
Pantelides' algorithm reports that the index is zero whereas it is indeed two.
This is because the method cannot detect the singularity of the coefficient matrix $\begin{pmatrix} 1 & 1 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 1\end{pmatrix}$ of $\begin{pmatrix} \dot{x}_1 \\ \dot{x}_2 \\ \dot{x}_3 \end{pmatrix}$.
As this toy example shows, Pantelides' method, which ignores numerical information, may fail on some DAEs due to numerical cancellations.
This kind of failure can also occur in other methods to reduce the index or to analyze DAEs such as a structural algorithm of Unger et al.~\cite{Unger1995} and the $\Sigma$-method of Pryce~\cite{Pryce2001}.

Some index reduction algorithms address this problem.
One example is the $\sigma \nu$-method by Chowdhry et al.~\cite{Chowdhry2004}, which is based on the algorithm by Unger et al.~\cite{Unger1995}.
The method performs the Gaussian elimination on the Jacobian matrix $\partial F / \partial \dot{x}$ under the assumption that nonlinear or time-varying terms do not cancel out.
Tan et al.~\cite{Tan2016} presented two methods called LC-method and ES-method, which improve the $\Sigma$-method of Pryce.
Their methods identify numeric or symbolic cancellations and modify the DAE if necessary.

Pantelides' method~\cite{Pantelides1988} and the $\Sigma$-method~\cite{Pryce2001} discard numerical information, which sometimes leads to a failure of the methods.
In dynamical systems, specific numbers in structural equations such as the conservation laws should be treated numerically, while we can deal with physical characteristic values as nonzero parameters without reference to their nominal values.
For a faithful model of a dynamical system, it is natural to distinguish accurate and inaccurate numbers.
This led Murota--Iri~\cite{Murota1985} to introduce the notion of a \emph{mixed matrix}, which is a matrix consisting of the following two kinds of entries:
\begin{description}
  \item[Accurate Constants\normalfont{,}]
    which represent precise values such as coefficients of conservation laws.
    We assume that arithmetic operations of these constants can be performed in constant time.

  \item[Independent Parameters\normalfont{,}]
    which are algebraically independent over the field of accurate constants.
    These parameters often represent physical quantities such as masses, lengths or electric resistances since their nominal values are inaccurate by measurement noises and other errors.
    These parameters should be treated combinatorially without reference to their nominal values.
\end{description}

For example, consider an electric network consisting of voltage sources, resistances and wires connecting them.
A linear equation representing the circuit has two kinds of coefficients: the exact `$\pm 1$'s coming from Kirchhoff's law, and the resistance values coming from Ohm's law.
Since the nominal values of resistances are supposed to be inaccurate, it is natural to model the system by a linear equation with a mixed matrix, where constants and parameters represent the exact `$\pm 1$'s and the resistances, respectively.
See an example in Section~\ref{sec:example_2} for modeling of an RLC circuit with a mixed matrix.

Mathematical aspects of mixed matrices are studied by Murota~\cite{Murota1993}.
If all nonzero entries of a matrix are independent parameters, then its rank is equal to the maximum size of a matching in an associated bipartite graph.
For a mixed matrix, the rank computation corresponds to solving an independent matching problem on matroids, which is a generalization of the maximum matching problem on bipartite graphs.
An efficient algorithm based on matroid theory is provided for the rank computation of mixed matrices; see \cite{Murota2000} for detail.

In this paper, we provide an index reduction algorithm for a linear DAE
\begin{align} \label{def:linear_dae}
  \sum_{k=0}^l A_k \app{x^{(k)}}{t} = \app{f}{t}
\end{align}
with $n \times n$ mixed matrices $A_0, A_1, \ldots, A_l$ and a sufficiently smooth function $\funcdoms{f}{\setR}{\setR^n}$.

Modeling dynamical systems by mixed matrices is advantageous for the index reduction as follows.
First, we can resolve numerical cancellations in a mixed matrix by modifying only accurate constants since the cancellations do not occur involving independent parameters.
Thus our algorithm is expected to preserve the sparsity of the matrix compared to other index reduction algorithms that transform the entire matrix.
Next, since accurate constants arising from typical dynamical systems are integers or rational numbers, we can avoid arithmetic operations and comparisons of floating-point numbers.
This fact makes the algorithm stable.
Finally, if we want to simulate a dynamical system on many different settings of physical quantities, we can reuse the resulting low-index DAE as long as the values of physical quantities do not unluckily cancel out.

Our approach for the index reduction is to transform a DAE into an equivalent DAE to which the MS-algorithm is applicable, based on the \emph{combinatorial relaxation} method, introduced by Murota~\cite{Murota1990, Murota1995} as a framework to solve a linear algebraic problem by iteratively relaxing it into an efficiently solvable combinatorial optimization problem.
Our algorithm does not rely on symbolic manipulations, while it uses numerical calculations and fast combinatorial optimization algorithms on graphs and matroids.
For an $l$-th order DAE, our algorithm is proved to run in $\Order{l^2 n^{\omega + 2}}$ time, where $\omega$ is the matrix multiplication exponent, i.e., the number of arithmetical operations needed to multiply two $n \times n$ matrices is $\Order{n^\omega}$.
The current best known value of $\omega$ is almost 2.3728639 due to~\cite{Gall2014}.
In practice, however, one can adopt $\omega = 3$.
Our algorithm is expected to run much faster in most cases because
the algorithm terminates without modifying a DAE unless it has unlucky numerical cancellations.

In addition, we give an improved algorithm for DAEs whose coefficients are \emph{dimensionally consistent}.
The dimensional consistency, which is introduced by Murota~\cite{Murota1985b}, is a mathematical assumption on mixed matrices reflecting the principle of dimensional homogeneity in physical systems, and DAEs arising from dynamical systems naturally ensure this assumption.
We show that the improved algorithm retains the dimensional consistency, and
that the running time is $\Order{l n^4 \log n}$.
Since the order $l$ of typical DAEs is not so large, we claim that the algorithm does not significantly increase the computational cost compared to the MS-algorithm, which costs $\Order{n^4 + ln^2}$ time by the straightforward matrix multiplication.

Furthermore, though our index reduction algorithm is designed for linear DAEs, the algorithm can be applied to nonlinear DAEs by regarding nonlinear terms as independent parameters.
While the similar approach is adopted in $\sigma \nu$-method of Chowdhry et al.~\cite{Chowdhry2004}, our method is expected to be applicable to a larger class of nonlinear DAEs than the $\sigma \nu$-method because our method does not transform DAEs involving nonlinear terms.

We describe the relation between the proposed algorithm and related index reduction algorithms.
If all nonzero entries of $A(s)$ are independent parameters, our algorithm just passes a given DAE to the MS-method as we described above.
In contrast, if $A(s)$ has no independent parameter, then our algorithm coincides with the LC-method by Tan et al.~\cite{Tan2016}.
In addition, our algorithm is similar to the $\sigma \nu$-method~\cite{Chowdhry2004} in the sense that both methods treat matrices having accurate constants and independent parameters, yet their approaches are quite different; the $\sigma \nu$-method is based on the Gaussian elimination approach by Gear~\cite{Gear1988}, whereas our algorithm relies on the dummy variable approach by Mattsson--S\"{o}derlind~\cite{Mattsson1993}.

A recent work~\cite{Iwata2017} has proposed an index reduction algorithm which is proved to work for any instances of first order linear DAEs with constant coefficients.
The algorithm directly reduces the index of a given DAE by row operations, whereas our algorithm only resolves numerical cancellations in a DAE and eventually relies on the MS-algorithm for the actual index reduction process.
Thus our algorithm is expected to preserve the sparsity of DAEs compared to the algorithm in~\cite{Iwata2017}.

This paper is organized as follows.
Section~\ref{sec:index_reduction_for_linear_DAEs} reviews the previous index computation and reduction algorithms for linear DAEs with constant coefficients, including the MS-algorithm and combinatorial relaxation algorithms.
Section~\ref{sec:mixed_matrices} explains mixed matrices and their rank identities.
Section~\ref{sec:algorithm} describes the proposed algorithm.
Section~\ref{sec:exploiting_dimensional_consistency} improves our algorithm under the assumption of the dimensional consistency.
Section~\ref{sec:example} illustrates two examples.
Section~\ref{sec:nonlinear_dae} discusses an application to nonlinear DAEs.
Finally, Section~\ref{sec:conclusion} concludes this paper.

\section{Index Reduction for Linear DAEs}
\label{sec:index_reduction_for_linear_DAEs}

\subsection{Index of Linear DAEs}

A linear DAE with constant coefficients is the following DAE:
\begin{align} \label{def:linear_dae2}
  \sum_{k=0}^l A_k \app{x^{(k)}}{t}= \app{f}{t},
\end{align}
where $A_0, A_1, \ldots, A_l$ are $n \times n$ matrices and $\funcdoms{f}{\setR}{\setR^n}$ is a sufficiently smooth function.
We assume that $f$ is Laplace transformable for simplicity, though this assumption is not essential.
By the Laplace transformation, the DAE~\eqref{def:linear_dae2} is transformed into
\begin{align} \label{eq:laplaced_dae2}
  A(s) \tilde{x}(s) = \tilde{f}(s) + \sum_{k=0}^l \sum_{i=1}^k s^{k-i} A_k x^{(i-1)}(0),
\end{align}
where $\tilde{x}(s)$ and $\tilde{f}(s)$ are the Laplace transforms of $x(t)$ and $f(t)$, respectively, and $A(s) = \sum_{k=0}^l s^k A_k$.
We henceforth denote the right-hand side of~\eqref{eq:laplaced_dae2} by $\hat{f}(s)$.
The matrix $A(s)$ is a matrix whose entries are polynomials, called a \emph{polynomial matrix}.
We say that $A(s)$ is \emph{nonsingular} if its determinant is not identically zero.

An initial value $\prn[\big]{x_0, x^{(1)}_0, \ldots, x^{(l-1)}_0} \in \setR^n \times \cdots \times \setR^n$ is said to be \emph{consistent} if there exists at least one solution of~\eqref{def:linear_dae2} satisfying
\begin{align} \label{cond:initial_value}
  x(0) = x_0, \, \dot{x}(0) = x^{(1)}_0, \, \ldots, \, x^{(l-1)}(0) = x^{(l-1)}_0.
\end{align}
We say that the DAE~\eqref{def:linear_dae2} is \emph{solvable} if there exists a unique solution of~\eqref{def:linear_dae2} satisfying the initial value condition~\eqref{cond:initial_value} for an arbitrary consistent point.
The solvability of~\eqref{def:linear_dae2} is characterized by $A(s)$ as follows.

\begin{theorem}[{\cite{Brenan1996, Shi2004}}]
  A linear DAE~\eqref{def:linear_dae2} is solvable if and only if the associated polynomial matrix $A(s)$ is nonsingular.
\end{theorem}

See~\cite[Theorem~2.3.1]{Brenan1996} for $l=1$ and~\cite[Theorems~2.22--23]{Shi2004} for $l \geq 2$.
In this paper, we focus on solvable DAEs~\eqref{def:linear_dae2}.
In addition, we abuse the term ``DAE'' and also refer to the equation~\eqref{eq:laplaced_dae2} as a DAE.

The differentiation index of a first order linear DAE~\eqref{def:linear_dae2} with $A(s) = A_0 + sA_1$ is known to equal
\begin{align} \label{def:nu_A}
  \nu(A) = \delta_{n-1}(A) - \degdet{A} + 1.
\end{align}
Here, $\delta_k(A)$ denotes the maximum degree of the determinant of a submatrix in $A(s)$ of size $k$, i.e.,
\begin{align*}
  \delta_k(A) = \max \set{\deg \det A(s)[I, J]}[\card{I} = \card{J} = k],
\end{align*}
where $A(s)[I, J]$ is the submatrix in $A(s)$ with row set $I$ and column set $J$, and $\deg p(s)$ designates the degree of a polynomial $p(s)$ in $s$.
In particular, $\degdet{A}$ is the degree of the determinant of $A(s)$, and  $\delta_{n-1}(A)$ is the maximum degree of a cofactor of $A(s)$.
For a DAE~\eqref{def:linear_dae2} with $l \geq 2$, its index is defined to be that of a first order DAE obtained by replacing higher order derivatives with new variables~\cite{Tan2016}.

\subsection{Assignment Problem}

In analysis of DAEs, Pryce~\cite{Pryce2001} introduced an assignment problem as a reinterpretation of Pantelides' algorithm~\cite{Pantelides1988} as follows.

Consider a linear DAE~\eqref{eq:laplaced_dae2} with $n \times n$ nonsingular polynomial matrix $\app{A}{s}$ with row set $R$ and column set $C$.
We denote the $\prn{i, j}$ entry of $A(s)$ by $A_{i,j}(s)$.
Let $\app{G}{A}$ denote the bipartite graph with vertex set $R \cup C$ and edge set $\app{E}{A} = \set{\prn{i, j} \in R \times C}[A_{i,j}(s) \neq 0]$.
An edge subset $M \subseteq \app{E}{A}$ is called a \emph{matching} if the ends of edges in $M$ are disjoint.
Since $A(s)$ is nonsingular, $G(A)$ has a matching of size $n$, called a \emph{perfect matching}.
We set the weight $c_e$ of an edge $e = \prn{i, j} \in \app{E}{A}$ by $c_e = c_{i,j} = \deg \app{A_{i,j}}{s}$.

The assignment problem on $G(A)$ is the following problem $\P{A}$:
\Maximize[name=$\P{A}$]{
  \sum_{e \in M} c_e
}{
  \text{$M \subseteq \app{E}{A}$ is a perfect matching on $\app{G}{A}$.}
}
The dual problem $\D{A}$ of $\P{A}$ is expressed as follows:
\Minimize[name=$\D{A}$]{
  \sum_{j \in C} q_j - \sum_{i \in R} p_i
}{
  q_j - p_i \geq c_{i, j} & \prn{\prn{i, j} \in \app{E}{A}}, \\
  p_i \in \setZ           & \prn{i \in R}, \\
  q_j \in \setZ           & \prn{j \in C}.
}

We denote the optimal value of the problem $\P{A}$ (and $\D{A}$) by $\tdegdet{A}$.
Recall that $\degdet{A}$ denotes $\deg \det \app{A}{s}$.
It is well known that $\degdet{A} \leq \tdegdet{A}$ holds, and the equality is attained if and only if the coefficient of $s^{\tdegdet{A}}$ in $\det \app{A}{s}$ does not vanish.
In this sense, $\tdegdet{A}$ serves as a combinatorial upper bound on $\degdet{A}$.
We call $\app{A}{s}$ \emph{upper-tight} if $\degdet{A} = \tdegdet{A}$ holds.

For a dual feasible solution $\prn{p, q}$, a \emph{tight coefficient matrix} $A^\#$ of $\app{A}{s}$ is defined by
\begin{align*}
  A_{i,j}^\# \defeq \text{the coefficient of $s^{q_j - p_i}$ in $A_{i,j}(s)$}
\end{align*}
for each $i \in R$ and $j \in C$.
Note that $A^\#$ changes depending on $\prn{p, q}$.
This matrix is called a ``$\Sigma$-Jacobian'' by Pryce~\cite{Pryce2001}; the name ``tight coefficient matrix'' is due to Murota~\cite{Murota1995b}.

\subsection{Computing the Index via Combinatorial Relaxation}
\label{sec:computing_the_index_via_the_combonatorial_relaxation}

The tight coefficient matrix plays an important role in the combinatorial relaxation algorithm of Murota~\cite{Murota1995} to compute $\delta_n(A)$ for a polynomial matrix $A(s)$ through the following lemma.

\begin{lemma}[{\cite[Proposition 6.2]{Murota1990}}] \label{lem:tightness1}
  Let $A(s)$ be an $n \times n$ nonsingular polynomial matrix and $A^\#$ the tight coefficient matrix of $\app{A}{s}$ with respect to an optimal solution of $\D{A}$.
  Then $\app{A}{s}$ is upper-tight if and only if $A^\#$ is nonsingular.
\end{lemma}

The combinatorial relaxation method consists of the following three phases.

\begin{enumerate}[label={Phase \arabic*.}]
  \item Compute a combinatorial upper bound $\hat{\delta}_{n}(A)$ of $\delta_{n}(A)$ by solving an assignment problem.
  \item Check whether $A(s)$ is upper-tight using Lemma~\ref{lem:tightness1}.
        If it is, return the estimation and halt.
  \item Modify $A(s)$ to improve $\hat{\delta}_{n}(A)$ by replacing $A(s)$ with $U(s)A(s)$, where $U(s)$ is a \emph{unimodular matrix}.
        Go back to Phase 2.
\end{enumerate}

Here, a unimodular matrix is a square polynomial matrix whose determinant is a nonzero constant.
The algorithm is designed so that $\tdegdet{A}$ decreases in each iteration, while unimodular transformations preserve $\degdet{A}$.
Thus, after a finite number of iterations, the algorithm terminates with $\tdegdet{A} = \degdet{A}$.

Subsequently, Murota~\cite{Murota1995b} applied the combinatorial relaxation approach to computing $\delta_k(A)$ for $k = 1, \ldots, n$.
In this algorithm, Phase~3 modifies $A(s)$ to $U(s)A(s)V(s)$, where $U(s)$ and $V(s)$ are \emph{biproper Laurent polynomial matrices}, i.e., entries are all polynomials in $1/s$ and the determinants are nonzero constants.
This type of transformation is known to preserve $\delta_{k}(A)$.
The values of $\delta_{n-1}(A)$ and $\delta_n(A)$ determine the index $\nu(A)$ by~\eqref{def:nu_A}.

\subsection{Mattsson--S\"{o}derlind's Index Reduction Algorithm}
\label{sec:ms_alg}

We now review Mattsson--S\"{o}derlind's index reduction algorithm applied to a linear DAE~\eqref{eq:laplaced_dae2} with $n \times n$ nonsingular polynomial matrix $A(s)$.
Let $\prn{p, q}$ be an optimal solution of $\D{A}$.
For $h \in \setZ$, we define
\begin{align*}
  R_h \defeq \set{i \in R}[p_i = h], &\quad R_{\geq h} \defeq \set{i \in R}[p_i \geq h], \\
  C_h \defeq \set{j \in C}[q_j = h], &\quad C_{\geq h} \defeq \set{j \in C}[q_j \geq h].
\end{align*}
The MS-algorithm applied to the DAE~\eqref{eq:laplaced_dae2} is outlined as follows.
The following description is a rewritten version especially for linear DAEs, though the original MS-algorithm is designed for nonlinear DAEs~\cite[Section 3.1]{Mattsson1993}.

\paragraph{Mattsson--S\"{o}derlind's Index Reduction Algorithm}
\begin{enumerate}[label={Step \arabic*.}]
  \item Compute an optimal solution $\prn{p, q}$ of $\D{A}$ satisfying $p_i, q_j \geq 0$ for $i \in R$ and $j \in C$.
        Let $A^\#$ denote the tight coefficient matrix of $A(s)$ with respect to $\prn{p, q}$.
        If $A^\#$ is singular, then the algorithm terminates in failure.
  \item For each $h = 0, \ldots, \eta + 1 \,\, \prn[\big]{\eta \defeq \displaystyle \max_{i \in R} p_i}$, obtain $J_h \subseteq C_{\geq h}$ such that $A^\#[R_{\geq h}, J_h]$ is nonsingular and
        %
        \begin{align*}
          C = J_0 \supseteq J_1 \supseteq J_2 \supseteq \cdots \supseteq J_{\eta} \supseteq J_{\eta + 1} = \varnothing.
        \end{align*}
  \item For each $j \in C$, let $k_j$ be an integer such that $j \in J_{k_j}$ and $j \notin J_{k_j + 1}$.
        Introduce $k_j$ dummy variables $z^{[q_j]}_j, z^{[q_j-1]}_j, \ldots, z^{[q_j-k_j+1]}_j$ corresponding to $s^{q_j}\tilde{x}_j, \, s^{q_j-1}\tilde{x}_j, \ldots, s^{q_j-k_j+1}\tilde{x}_j$, respectively.
  \item For each $i \in R$, return the 0-th, 1-st, ..., $p_i$-th order derivatives of the $i$-th equation.
        Replace variables with the corresponding dummy variables.
\end{enumerate}

The validity of the MS-algorithm is established as follows.

\begin{proposition}[{\cite[Section 3.2]{Mattsson1993}}] \label{prop:ms-alg}
  Let $A(s)$ be a polynomial matrix in the DAE~\eqref{eq:laplaced_dae2} and $A^\#$ the tight coefficient matrix of $A(s)$ with respect to an optimal solution of $\D{A}$.
  If $A^\#$ is nonsingular, then the MS-algorithm correctly returns an equivalent DAE with index at most one.
\end{proposition}

From Lemma~\ref{lem:tightness1}, the condition in Proposition~\ref{prop:ms-alg} is equivalent to the upper-tightness of $A(s)$ as follows.

\begin{corollary} \label{cor:MS_equivalent_validity_condition}
  Let $A(s)$ be a polynomial matrix in the DAE~\eqref{eq:laplaced_dae2}.
  If $A(s)$ is upper-tight, then the MS-algorithm correctly returns an equivalent DAE with index at most one.
\end{corollary}

The description above is still valid for a nonlinear DAE~\eqref{def:dae} by redefining $A(s)$ as
\begin{align*}
  \text{the coefficient of $s^k$ in $A_{i,j}(s)$} \defeq \text{the partial derivative of the $i$-th equation with respect to $x^{(k)}_j$}
\end{align*}
for each $i = 1, \ldots, n$, $j = 1, \ldots, n$ and $k = 0, \ldots, l$.
Then the nonsingularity of $A^\#$ essentially comes from the requirement of the implicit function theorem, which is used to convert the DAE into an ODE by solving the DAE for the highest order derivatives.

\subsection{Index Reduction via Combinatorial Relaxation}
\label{sec:index_reduction_algorithm}

For a linear DAE~\eqref{eq:laplaced_dae2} that does not satisfy the validity condition of the MS-algorithm, we need to modify it to apply the MS-algorithm.
Here, the modification of DAEs must preserve the sets of their solutions.
We can use unimodular transformations in the form of
\begin{align*}
  U(s) A(s) \tilde{x}(s) = U(s)\hat{f}(s),
\end{align*}
where $U(s)$ is a unimodular matrix.
Since unimodular transformations correspond to the operations of adding an equation in the DAE or its (higher order) derivative to another equation, the DAEs before and after the transformation have the same solution set.

Murota's combinatorial relaxation algorithm~\cite{Murota1995} for computing $\delta_n(A)$ described in Section~\ref{sec:computing_the_index_via_the_combonatorial_relaxation} modifies a given polynomial matrix $A(s)$ into an upper-tight polynomial matrix $\bar{A}(s) = U(s)A(s)$ using some unimodular matrix $U(s)$.
Then from Corollary~\ref{cor:MS_equivalent_validity_condition}, the matrix $\bar{A}(s)$ satisfies the validity condition of the MS-algorithm.
Therefore, we can use Murota's algorithm as an index reduction algorithm by combining with the MS-algorithm.
Note that this modification may change $\delta_{n-1}(A)$, and hence $\nu(A)$.
We also remark that this index reduction algorithm indeed coincides with the LC-method of Tan et al.~\cite{Tan2016} applied to the linear DAEs with constant coefficients.

\section{DAEs with Mixed Matrices}
\label{sec:mixed_matrices}

The algorithms explained in Section~\ref{sec:index_reduction_for_linear_DAEs} work under the assumption that one can exactly perform arithmetic operations of inaccurate physical quantities, which naturally arise from dynamical systems.
In order to distinguish inaccurate quantities from accurate numbers, we deal with a DAE with mixed matrices.

\subsection{Mixed Matrices and Mixed Polynomial Matrices}
Let $\setF$ be a field and $\setK$ a subfield of $\setF$.
A typical setting in the context of DAEs is $\setF = \setR$ and $\setK = \setQ$.
A matrix $T$ over $\setF$ is said to be \emph{generic} if the set  of nonzero entries of $T$ is algebraically independent over $\setK$.
A \emph{mixed matrix} with respect to $\prn{\setK, \setF}$ is a matrix in the form of $Q+T$, where $Q$ is a matrix over $\setK$ and $T$ is a generic matrix.
A \emph{layered mixed matrix} (or \emph{LM-matrix}) is a mixed matrix such that nonzero rows of $Q$ and $T$ are disjoint.
An LM-matrix $A$ can be expressed as $A = \binom{Q}{T}$.

A polynomial matrix $A(s) = \sum_{k=0}^l s^k A_k$ is called a \emph{mixed polynomial matrix} if it is expressed as $A_k = Q_k + T_k$ with $Q_k$ and $T_k$ that satisfy the following conditions:
\begin{enumerate}[label={(MP-Q)}]
  \item[(MP-Q)] Each $Q_k$ $\prn{k = 0, \ldots, l}$ is a matrix over $\setK$.
  \item[(MP-T)] The set of nonzero entries of $T_0, \ldots, T_l$ is algebraically independent over $\setK$.
\end{enumerate}
A \emph{layered mixed polynomial matrix} (or \emph{LM-polynomial matrix}) is a mixed polynomial matrix such that nonzero rows of $Q(s) = \sum_{k=0}^l s^k Q_k$ and $T(s) = \sum_{k=0}^l s^k T_k$ are disjoint.
An LM-polynomial matrix is expressed as $A(s) = \binom{Q(s)}{T(s)}$.

\subsection{Rank of LM-matrices}
For a matrix $A$, we denote the row set and column set by $\row{A}$ and $\col{A}$, respectively.
Consider the associated bipartite graph $\app{G}{A} = \prn{R, C; \app{E}{A}}$, where $R = \row{A}$ and $C = \col{A}$.
The \emph{term-rank} of $A$ is the maximum size of a matching in $G(A)$, and is denoted by $\trank A$.
It is well known that $\rank A \leq \trank A$ holds.
The equality is attained if and only if $A$ has a submatrix of size $\trank A$ with nonzero determinant.
This is analogous to the relation between $\delta_n$ and $\hat{\delta}_n$ for a polynomial matrix.

Let $A = \binom{Q}{T}$ be an LM-matrix.
If $A$ has no accurate constants, i.e., $A$ is a generic matrix $T$, it holds that $\rank T = \trank T$ from the independence of nonzero entries.
From this equality, we can compute $\rank T$ by solving a maximum matching problem on the associated bipartite graph $G(T)$.
For general LM-matrices, the following holds from the generalized Laplace expansion.
\begin{proposition}[{\cite[Theorem 3.1]{Murota1987}}]
  For an LM-matrix $A = \binom{Q}{T}$ with $R_Q = \row{Q}$, $R_T = \row{T}$ and $C = \col{A}$, the following rank identity holds:
  \begin{align}
    \label{eq:rank_identity_max} \rank A &= \max \set{\rank Q[R_Q, J] + \trank T[R_T, C \setminus J]}[J \subseteq C].
  \end{align}
\end{proposition}

The problem of maximizing the right-hand side of~\eqref{eq:rank_identity_max} can be reduced to an \emph{independent matching problem} on a matroid; see~\cite[Section~4.2]{Murota2000} for detail.
The following identity is obtained from the duality of the independent matching problem.
\begin{proposition}[{\cite[Theorem 3.1]{Murota1987}}]
  For an LM-matrix $A = \binom{Q}{T}$ with $R_Q = \row{Q}$, $R_T = \row{T}$ and $C = \col{A}$, the following rank identity holds:
  \begin{align}
    \label{eq:rank_identity_min} \rank A &= \min \set{\rank Q[R_Q, J] + \trank T[R_T, J] + \card{C \setminus J}}[J \subseteq C].
  \end{align}
\end{proposition}

Similarly, we give the following term-rank identity for LM-matrices, which will be used later in the proof of Lemma~\ref{lem:tilde_A_optimality}.

\begin{proposition} \label{prop:trank_identity_min}
  For an LM-matrix $A = \binom{Q}{T}$ with $R_Q = \row{Q}$, $R_T = \row{T}$ and $C = \col{A}$, the following term-rank identity holds:
  \begin{align*}
    \trank A &= \min \set{\trank Q[R_Q, J] + \trank T[R_T, J] + \card{C \setminus J}}[J \subseteq C].
  \end{align*}
\end{proposition}
\begin{proof}
  This immediately follows from the well-known rank formula of a union matroid~\cite{Edmonds1968} and the fact that the union of transversal matroids is also a transversal matroid~\cite[Corollary~11.3.8]{Oxley2011}.
\end{proof}

\subsection{Combinatorial Relaxation Algorithm for Mixed Polynomial Matrices}
\label{ref:combinatorial_relaxation_algorithm_for_mixed_polynomial_matrices}

Murota~\cite{Murota1998} described the first algorithm to compute $\delta_k(A)$ for a mixed polynomial matrix $A(s)$ through a reduction to a \emph{valuated independent assignment problem}.
A combinatorial relaxation algorithm for the computation of $\delta_k(A)$ is given in~\cite{Iwata2013}.
This algorithm first converts a mixed polynomial matrix into an LM-polynomial matrix $A(s) = \binom{Q(s)}{T(s)}$, and modifies $A(s)$ to
\begin{align*}
  \bar{A}(s) =
  \begin{pmatrix} U_Q(s) & O \\ O & I \end{pmatrix}
  \binom{Q(s)}{T(s)},
\end{align*}
where $I$ is an identity matrix of appropriate size and $U_Q(s)$ is a \emph{Laurent polynomial matrix}.
With the use of~\eqref{def:nu_A}, we can obtain the index $\nu(A)$ by computing $\delta_n(A)$ and $\delta_{n-1}(A)$.

In order to devise an index reduction algorithm for DAEs with mixed matrices, we need to make use of unimodular transformations instead of Laurent polynomial transformations, as explained in Section~\ref{sec:index_reduction_algorithm}.

\section[%
  Combinatorial Relaxation Algorithm for Index Reduction with Mixed Polynomial Matrices%
]{%
  Combinatorial Relaxation Algorithm for Index Reduction\\%
  with Mixed Polynomial Matrices%
}
\label{sec:algorithm}

This section presents our index reduction algorithm for a DAE
\begin{align} \label{eq:laplaced_mixed_dae}
  A(s) \tilde{x}(s) = \hat{f}(s)
\end{align}
with a nonsingular mixed polynomial matrix $A(s)$, which is the Laplace transform of the DAE~\eqref{def:linear_dae}.
From Corollary~\ref{cor:MS_equivalent_validity_condition}, our goal is to find a unimodular matrix $U(s)$ such that $\bar{A}(s) = U(s) A(s)$ is upper-tight.
Then applying the MS-algorithm to the DAE $U(s) A(s) \tilde{x}(s) = U(s) \hat{f}(s)$, we obtain a resultant low-index DAE.

We cannot perform row operations on $A(s)$ involving rows containing independent parameters.
Our first step is to convert a given DAE~\eqref{eq:laplaced_mixed_dae} into another DAE whose coefficient matrix $A(s)$ is an LM-polynomial matrix expressed as $A(s) = \binom{Q(s)}{T(s)}$.
Then we can transform $A(s)$ to
\begin{align} \label{eq:update_A}
  \bar{A}(s)
  =
  \begin{pmatrix} U_Q(s) & O \\ O & I \end{pmatrix}
  \binom{Q(s)}{T(s)},
\end{align}
where $U_Q(s)$ is a unimodular matrix.
Note that we are allowed to perform row operations only on $Q(s)$ even for an LM-polynomial matrix $A(s) = \binom{Q(s)}{T(s)}$, and thus we cannot always reduce the index to one only by row operations on $Q(s)$.
We describe this conversion process from mixed polynomial matrices into LM-polynomial matrices in Section~\ref{sec:mixed_to_LM}.

After the conversion, we find a unimodular matrix $U_Q(s)$ in~\eqref{eq:update_A} such that $A(s)$ is upper-tight based on the combinatorial relaxation approach.
The outline of our algorithm is shown as follows.

\paragraph{Algorithm for Tightness}
\begin{enumerate}[label={Phase \arabic*.}]
  \item Construct an optimal solution $\prn{p, q}$ of $\D{A}$.
  \item If the tight coefficient matrix $A^\#$ with respect to $\prn{p, q}$ is nonsingular, then return $A(s)$ and halt.
  \item Modify $A(s)$ into $\bar{A}(s)$ such that $\tdegdet{\bar{A}} \leq \tdegdet{A} - 1$ and $\degdet{A} = \degdet{\bar{A}}$.
        Update $\prn{p, q}$ to an optimal solution of $\D{\bar{A}}$, and go back to Phase~2.
\end{enumerate}
Section~\ref{sec:construct_dual_opt} describes an algorithm to find $\prn{p, q}$ in Phase~1.
The condition in Phase~2, which is equivalent to the upper-tightness of $A(s)$ by Lemma~\ref{lem:tightness1}, can be checked by solving an independent matching problem~\cite{Murota1987}.
The matrix modification and an update procedure of $\prn{p, q}$ in Phase 3 is explained in Sections~\ref{sec:matrix_modification} and~\ref{sec:dual_updates}, respectively.
In Section~\ref{sec:complexity_analysis}, we analyze the time complexity of our algorithm.

\subsection{Reduction to LM-polynomial Matrices}
\label{sec:mixed_to_LM}

We first convert the DAE~\eqref{eq:laplaced_mixed_dae} with a mixed polynomial coefficient matrix $A(s) = Q(s) + T(s)$ into the following augmented DAE
\begin{align} \label{eq:LM_dae}
  \begin{pmatrix}
     I &  Q(s) \\
    -D & DT(s)
  \end{pmatrix}
  \begin{pmatrix}
    \tilde{y}(s) \\
    \tilde{z}(s)
  \end{pmatrix}
  =
  \begin{pmatrix}
    \app{\hat{f}}{s} \\
    0
  \end{pmatrix},
\end{align}
where $D$ is a diagonal matrix whose diagonal entries are independent parameters $\tau_1, \ldots, \tau_n$.
Note that the coefficient matrix of the augmented DAE~\eqref{eq:LM_dae} is an LM-polynomial matrix as the set of nonzero coefficients of entries in $-D$ and $DT(s)$ is algebraically independent over $\setK$.

\begin{proposition}
  Let $\binom{\tilde{y}(s)}{\tilde{z}(s)}$ be a solution of the DAE~\eqref{eq:LM_dae}.
  Then $\tilde{z}(s)$ is a solution of the DAE~\eqref{eq:laplaced_mixed_dae}.
\end{proposition}
\begin{proof}
  By left-multiplying both sides of~\eqref{eq:LM_dae} by a nonsingular constant matrix $\begin{pmatrix} I & O \\ I & D^{-1} \end{pmatrix}$, we obtain
  %
  \begin{align*}
    \begin{pmatrix}
      I & Q(s) \\
      O & A(s)
    \end{pmatrix}
    \begin{pmatrix}
      \tilde{y}(s) \\
      \tilde{z}(s)
    \end{pmatrix}
    =
    \begin{pmatrix}
      \app{\hat{f}}{s} \\
      \app{\hat{f}}{s}
    \end{pmatrix},
  \end{align*}
  %
  where $O$ is a zero matrix.
  Thus it holds $A(s) \tilde{z}(s) = \hat{f}(s)$, which implies that $\tilde{z}(s)$ is a solution of the DAE~\eqref{eq:laplaced_mixed_dae}.
\end{proof}

\begin{example}
  Consider an index-3 DAE
  %
  \begin{align} \label{eq:dae_example_2}
    \begin{pmatrix}
      s & s & \alpha_1  \\
      s & s & 0 \\
      0 & \alpha_2 & \alpha_3s^2 -2s + \alpha_4
    \end{pmatrix}
    \begin{pmatrix}
      \tilde{x}_1(s) \\
      \tilde{x}_2(s) \\
      \tilde{x}_3(s)
    \end{pmatrix}
    =
    \begin{pmatrix}
      \hat{f}_1(s) \\
      \hat{f}_2(s) \\
      \hat{f}_3(s)
    \end{pmatrix},
  \end{align}
  %
  where $\alpha_1, \ldots, \alpha_4$ are independent parameters.
  Following~\eqref{eq:LM_dae}, we convert this DAE into
  %
  \begin{align} \label{eq:eae_example_2_LM}
    \prn{\begin{array}{ccc|ccc}
      1 & 0 & 0 & s & s & 0 \\
      0 & 1 & 0 & s & s & 0 \\
      0 & 0 & 1 & 0 & 0 & -2s \\
      \hline
      -\tau_1 &    0 &    0 & 0 & 0 & \alpha'_1 \\
      0    & -\tau_2 &    0 & 0 & 0 & 0 \\
      0    &    0 & -\tau_3 & 0 & \alpha'_2 & \alpha'_3 s^2 + \alpha'_4 \\
    \end{array}}
    \begin{pmatrix}
      \tilde{y}_1(s) \\
      \tilde{y}_2(s) \\
      \tilde{y}_3(s) \\
      \tilde{z}_1(s) \\
      \tilde{z}_2(s) \\
      \tilde{z}_3(s) \\
    \end{pmatrix}
    =
    \begin{pmatrix}
      \hat{f}_1(s) \\
      \hat{f}_2(s) \\
      \hat{f}_3(s) \\
      0 \\
      0 \\
      0 \\
    \end{pmatrix},
  \end{align}
  %
  where $\tau_1, \tau_2, \tau_3$ are independent parameters, $\alpha'_1 = \tau_1 \alpha_1$ and $\alpha'_i = \tau_3 \alpha_i$ for $i=2,3,4$.
  Then we can obtain a solution $\prn{\tilde{x}_1(s), \tilde{x}_2(s), \tilde{x}_3(s)}$ of~\eqref{eq:dae_example_2} by solving the augmented DAE~\eqref{eq:eae_example_2_LM}.
  While the index of~\eqref{eq:eae_example_2_LM} is also three, in general this conversion does not preserve the index of DAEs.
\end{example}

After the index reduction process, we need to fill independent parameters by real numbers to start numerical methods.
Indeed, we can substitute 1 for each diagonal entry $\tau_i$ of $D$, i.e., $D=I$.
To explain this fact, let
\begin{align} \label{def:B}
  B(s) = \begin{pmatrix} Q_1(s) & Q_2(s) \\ -D & DT(s) \end{pmatrix}
\end{align}
be the coefficient matrix of a DAE that our algorithm returns for the augmented DAE~\eqref{eq:LM_dae}, where $Q_1(s)$ and $Q_2(s)$ are some polynomial matrices.
By substituting the identity matrix to $D$, we obtain
\begin{align} \label{def:tilde_B}
  \bar{B}(s) = \begin{pmatrix} Q_1(s) & Q_2(s) \\ -I & T(s) \end{pmatrix}.
\end{align}
Though $\bar{B}(s)$ is no longer an LM-polynomial matrix, the following lemma guarantees the upper-tightness of $\bar{B}(s)$.
\begin{lemma} \label{lem:B}
  Let $Q_1(s)$, $Q_2(s)$ and $T(s)$ are polynomial matrices and $D$ a nonsingular diagonal matrix.
  Then $B(s)$ in~\eqref{def:B} is upper-tight if and only if $\bar{B}(s)$ in~\eqref{def:tilde_B} is upper-tight.
\end{lemma}
\begin{proof}
  Using $P = \begin{pmatrix} I & O \\ O & D^{-1} \end{pmatrix}$, we have $\bar{B}(s) = P B(s)$.
  Since $P$ is a nonsingular constant matrix, $\degdet{B} = \degdet{\bar{B}}$ holds.
  In addition, since $P$ is nonsingular, diagonal and constant, the row transformation by $P$ does not change the bipartite graph $G(B)$ and its edge weight $c_e$ associated with $B(s)$.
  This fact implies that $\tdegdet{B} = \tdegdet{\bar{B}}$.
  Thus the upper-tightness of $B(s)$ and $\bar{B}(s)$ are equivalent.
\end{proof}
From this lemma, we can ``forget'' the existence of $D$ in the augmented DAE~\eqref{eq:LM_dae}.
That is, to reduce the index of the DAE~\eqref{eq:laplaced_mixed_dae}, it suffices to apply our algorithm to a DAE
\begin{align*}
  \begin{pmatrix}
     I & Q(s) \\
    -I & T(s)
  \end{pmatrix}
  \begin{pmatrix}
    \tilde{y}(s) \\
    \tilde{z}(s)
  \end{pmatrix}
  =
  \begin{pmatrix}
    \app{\hat{f}}{s} \\
    0
  \end{pmatrix}
\end{align*}
as if the set of nonzero coefficients of entries in $\begin{pmatrix} -I & T(s) \end{pmatrix}$ were independent.

\subsection{Construction of Dual Optimal Solution}
\label{sec:construct_dual_opt}

For an LM-polynomial matrix $A(s) = \binom{Q(s)}{T(s)}$, an optimal solution $\prn{p, q}$ of $\D{A}$ is constructed as follows.
First, we obtain a maximum-weight perfect matching $M$ in $G(A)$ by the Hungarian method~\cite{Kuhn1955}.
Next, construct a residual graph $G_M = (W, E_M)$ with $W = R \cup C \cup \set{r}$ and $E_M = \app{E}{A} \cup \reor{M} \cup Z$, where $R = \row{A}$, $C = \col{A}$, $r$ is a new vertex, $\reor{M} = \set{\prn{j, i}}[\prn{i, j} \in M]$, and $Z = \set{\prn{r, i}}[i \in R]$.
The arc length $\funcdoms{\gamma}{E_M}{\setZ}$ of $G_M$ is defined by
\begin{align*}
  \app{\gamma}{i, j} = \begin{cases}
    -c_{i, j} & \prn{\prn{i, j} \in \app{E}{A}}, \\
     c_{j, i} & \prn{\prn{i, j} \in \reor{M}}, \\
     0        & \prn{\prn{i, j} \in Z}, \\
  \end{cases}
\end{align*}
for each $\prn{i, j} \in E_M$.
\begin{lemma}
  For a residual graph $G_M$ defined above, the following hold.
  %
  \begin{enumerate}[label={(\arabic*)}]
    \item All vertices are reachable from $r$ on $G_M$.
    \item There is no negative-weight directed cycle with respect to the arc length $\gamma$ in $G_M$.
  \end{enumerate}
\end{lemma}
\begin{proof}
  (1)
  Every vertex $i \in R$ is reachable from $r$ through an edge $\prn{r, i} \in Z$.
  In addition, since $\app{G}{A}$ has a perfect matching $M$, every vertex $j \in C$ is also reachable from $r$ via $i \in R$ through edges $\prn{r, i} \in Z$ and $\prn{i, j} \in M \subseteq \app{E}{A}$.

  (2)
  Suppose that $G_M$ has a simple directed cycle $\mathcal{C} \subseteq E_M$ of negative weight.
  By the construction of the residual graph $G_M$, we can express $\mathcal{C}$ as
  %
  \begin{align*}
    \mathcal{C} = \set{\prn{i_0, j_1}, \prn{j_1, i_1}, \prn{i_1, j_2}, \ldots, \prn{i_{L-1}, j_L}, \prn{j_L, i_L}},
  \end{align*}
  %
  where $L = \card{\mathcal{C}} / 2$, $i_0, \ldots, i_L \in R$, $j_1, \ldots, j_L \in C$ and $i_0 = i_L$.
  Notice that $\prn{i_{k-1}, j_k} \notin M$ and $\prn{j_k, i_k} \in M$ hold for $k = 1, \ldots, L$.
  Let
  %
  \begin{align*}
    M' \defeq \prn{M \setminus \set{\prn{j_1, i_1}, \ldots, \prn{j_L, i_L}}} \cup \set{\prn{i_0, j_1}, \ldots, \prn{i_{L-1}, j_L}}.
  \end{align*}
  %
  Now $M'$ is a perfect matching of $\app{G}{A}$, and since
  %
  \begin{align*}
    \sum_{e \in \mathcal{C}} \app{\gamma}{e} = -\sum_{k=1}^L c_{i_{k-1}, j_k} + \sum_{k=1}^L c_{j_k, i_k} < 0,
  \end{align*}
  %
  the weight of $M'$ is greater than the weight of $M$.
  This contradicts the assumption that $M$ is a maximum-weight perfect matching.
\end{proof}
This lemma guarantees the existence of a shortest path from $r$ to each vertex with respect to the arc length $\gamma$ in $G_M$.
Let $\app{d}{i, j}$ be the length of a shortest path from $i \in W$ to $j \in W$.
Using $d$, we define
\begin{align}
  p_i &\defeq -\app{d}{r, i} + \max_{i^* \in R} \app{d}{r, i^*}, \label{def:initial_p} \\
  q_j &\defeq -\app{d}{r, j} + \max_{i^* \in R} \app{d}{r, i^*}  \label{def:initial_q}
\end{align}
for each $i \in R$ and $j \in C$.

\begin{lemma} \label{lem:bound_pq}
  Let $\prn{p, q}$ be defined in~\eqref{def:initial_p} and~\eqref{def:initial_q}.
  Then $\prn{p, q}$ is an optimal solution of $\D{A}$ satisfying $0 \leq p_i \leq ln$ for each $i \in R$ and $0 \leq q_j \leq ln$ for each $j \in C$, where $n = \card{R} = \card{C}$ and $l$ is the maximum degree of an entry in $A(s)$.
\end{lemma}
\begin{proof}
  First, we prove that $\prn{p, q}$ is a feasible solution of $\D{A}$.
  By the definition of $\prn{p, q}$, every $p_i \, \prn{i \in R}$ and $q_j \, \prn{j \in C}$ are clearly integer.
  For each $\prn{i, j} \in \app{E}{A}$, it holds $\app{d}{r, j} \leq \app{d}{r, i} - c_{i, j}$.
  Thus
  %
  \begin{align*}
    q_j - p_i = -\app{d}{r, j} + \app{d}{r, i} \geq c_{i, j}
  \end{align*}
  %
  and this implies that $\prn{p, q}$ is a feasible solution of $\D{A}$.

  We second show the optimality of $\prn{p, q}$.
  For each $\prn{i, j} \in M$, since $\prn{i, j} \in E_M$ and $\prn{j, i} \in E_M$, we obtain
  %
  \begin{align*}
    q_j - p_i = -\app{d}{r, j} + \app{d}{r, i} = c_{i, j}.
  \end{align*}
  %
  Thus it holds that
  %
  \begin{align*}
    \sum_{j \in C} q_j - \sum_{i \in R} p_i = -\sum_{j \in C} \app{d}{r, j} + \sum_{i \in R} \app{d}{r, i} = \sum_{\prn{i, j} \in M} \prn{-\app{d}{r, j} + \app{d}{r, i}} = \sum_{\prn{i, j} \in M} c_{i,j}
  \end{align*}
  %
  which implies that $\prn{p, q}$ is optimal to $\D{A}$.

  Finally, we give the lower and upper bounds of each $p_i$ and $q_j$.
  The non-negativity of $p_i$ clearly holds from the definition of $p_i$.
  In addition, since $G(A)$ has a perfect matching, each $j \in C$ is incident to at least one vertex $i \in R$.
  Thus we obtain $q_j \geq p_i + c_{i,j} \geq 0$ by $p_i \geq 0$ and $c_{i, j} \geq 0$.
  Let $i^* \in R$ denote a vertex such that $\app{d}{r, i^*} \geq \app{d}{r, i}$ for each $i \in R$.
  Fix $j \in C$.
  Let $P_j \subseteq E_M$ and $P_{i^*} \subseteq E_M$ be shortest paths from $r$ to $j$ and $i^*$, respectively.
  Let $v \in W$ be the last common vertex in $P_j$ and $P_{i^*}$.
  Note that it holds
  %
  \begin{align*}
    q_j = -\app{d}{r, j} + \app{d}{r, i^*} = -\app{d}{v, j} + \app{d}{v, i^*}.
  \end{align*}
  %
  Let $Q_j \subseteq P_j$ and $Q_{i^*} \subseteq P_{i^*}$ be subpaths from $v$ to $j$ and $i^*$, respectively.
  From the definition of the arc lengths, $-\app{d}{v, j} \leq l \card{Q_j \cap \app{E}{A}}$ and $\app{d}{v, i^*} \leq l \card{Q_{i^*} \cap \reor{M}}$ hold.
  If $v = r$, we obtain
  %
  \begin{align*}
    q_j &\leq l\prn{\card{Q_j \cap \app{E}{A}} + \card{Q_{i^*} \cap \reor{M}}} \\
        &=    l\prn{\frac{\card{Q_j}}{2} + \frac{\card{Q_{i^*}} - 1}{2}} \\
        &< ln.
  \end{align*}
  %
  The last inequality holds since $\card{Q_j} + \card{Q_{i^*}} \leq \card[\big]{W \setminus \set{r}} = 2n$.
  Otherwise, it holds
  %
  \begin{align*}
    q_j &\leq l \prn{\card{Q_j \cap \app{E}{A}} + \card{Q_{i^*} \cap \reor{M}}} \\
        &\leq l \prn{\floor{\frac{\card{Q_j} + 1}{2}} + \floor{\frac{\card{Q_{i^*}} + 1}{2}}} \\
        &\leq ln,
  \end{align*}
  %
  where $\card{Q_j} + \card{Q_{i^*}} \leq \card[\big]{W \setminus \set{r, v}} = 2n - 1$ is used in the last inequality.
  Thus $q_j \leq ln$ holds for each $j \in C$.
  In addition, for each $i \in R$, we have $p_i \leq q_j - c_{i,j} \leq q_j \leq ln$, where $j \in C$ is incident to $i$ in $M$.
\end{proof}

\subsection{Matrix Modification}
\label{sec:matrix_modification}

Let $A(s) = \binom{Q(s)}{T(s)}$ be an $n \times n$ nonsingular LM-polynomial matrix that is not upper-tight.
Let $A^\# = \binom{Q^\#}{T^\#}$ be the tight coefficient matrix with respect to an optimal solution $\prn{p, q}$ of $\D{A}$.
Without loss of generality, we assume that $\row{Q} = R_Q = \set{1, \ldots, m_Q}$ and $p_1 \leq \cdots \leq p_{m_Q}$, where $m_Q = \card{R_Q}$.

Recall the rank identity~\eqref{eq:rank_identity_min}.
Let $J^* \subseteq C$ be a column subset that minimizes the right-hand side of the identity for $A^\#$, i.e., it holds
\begin{align} \label{eq:Jstar}
  \rank A^\# = \rank Q^\# [R_Q, J^*] + \trank T^\# [R_T, J^*] + \card{C \setminus J^*}.
\end{align}
Such $J^*$ is called a $\emph{minimizer}$ of~\eqref{eq:rank_identity_min}.
By the row transformation of $Q^\#$, we obtain a matrix $\bar{Q}^\# = UQ^\#$ such that
\begin{align} \label{eq:rank_Qsharp}
  \rank \bar{Q}^\#[R_Q, J^*] = \trank \bar{Q}^\#[R_Q, J^*].
\end{align}
In particular, this transformation can be accomplished only by operations of adding a scalar multiple of a row $i \in R_Q$ to another row $j \in R_Q$ with $p_i > p_j$.
Then the matrix $U$ is upper-triangular due to the order of rows in $R_Q$.
This is known as the \textit{backward elimination} of $\bar{Q}^\#[R_Q, J^*]$.
Consider
\begin{align} \label{def:U_tilde}
  U_Q(s) =
  \diag{s^{-p_1}, \ldots, s^{-p_{m_Q}}}
  U
  \diag{s^{p_1}, \ldots, s^{p_{m_Q}}},
\end{align}
where $\diag{a_1, \ldots, a_n}$ denotes a diagonal matrix with diagonal entries $a_1, \ldots, a_n$.
Note that each entry in $U_Q(s)$ is a polynomial because $U$ is upper-triangular.
In addition, since $\det U_Q(s) = \det U$ is a non-zero constant, $U_Q(s)$ is unimodular.

We define $D_p(s) = \diag{s^{p_1}, \ldots, s^{p_n}}$ and $D_q(s)=\diag{s^{q_1}, \ldots, s^{q_n}}$.
Using $U_Q(s)$, we update $A(s)$ to $\bar{A}(s)$ as in~\eqref{eq:update_A}:
\begin{align} \label{eq:update_A_2}
  \bar{A}(s) =
  \begin{pmatrix}
    U_Q(s) & O \\
    O      & I
  \end{pmatrix}
  A(s)
  =
  D_p^{-1}(s)
  \begin{pmatrix}
    U & O \\
    O & I
  \end{pmatrix}
  D_p(s)
  A(s).
\end{align}

To show that $\prn{p, q}$ is not an optimal solution of $\D{\bar{A}}$, we use the following lemma, which is given by Murota~\cite{Murota1990} as a combinatorial counterpart to Lemma~\ref{lem:tightness1}.

\begin{lemma}[{\cite[Proposition 6.2]{Murota1990}}] \label{lem:tightness2}
  Let $A(s)$ be an $n \times n$ nonsingular polynomial matrix and $A^\#$ the tight coefficient matrix of $\app{A}{s}$ with respect to a feasible solution $\prn{p, q}$ of $\D{A}$.
  Then $\prn{p, q}$ is optimal if and only if $\trank A^\# = n$.
\end{lemma}

\begin{lemma} \label{lem:tilde_A_optimality}
  Let $A(s) = \binom{Q(s)}{T(s)}$ be an $n \times n$ nonsingular LM-polynomial matrix that is not upper-tight, and $A^\# = \binom{Q^\#}{T^\#}$ the tight coefficient matrix with respect to an optimal solution $\prn{p, q}$ of $\D{A}$.
  Then for the LM-polynomial matrix $\bar{A}(s)$ defined in~\eqref{eq:update_A_2}, the value $\prn{p, q}$ is feasible on $\D{\bar{A}}$ but not optimal.
\end{lemma}
\begin{proof}
  Consider a rational function matrix
  %
  \begin{align} \label{def:F}
    H(s) = D_p(s) \bar{A}(s) D_q^{-1}(s).
  \end{align}
  %
  For each $i \in R$ and $j \in C$, it holds that $\deg H_{i,j}(s) = \bar{c}_{i, j} + p_i - q_j$, where $\bar{c}_{i,j} = \deg \bar{A}_{i,j}(s)$.
  By substituting~\eqref{eq:update_A_2} into~\eqref{def:F}, we obtain
  %
  \begin{align*}
    H(s) &= \begin{pmatrix}
      U & O \\
      O & I
    \end{pmatrix} D_p(s) A(s) D_q^{-1}(s)
    = \begin{pmatrix}
      U & O \\
      O & I
    \end{pmatrix} \prn{A^\# + A^\infty(s)},
  \end{align*}
  %
  where $A^\infty(s)$ is a matrix whose nonzero entries are polynomials in $s^{-1}$.
  Hence for each $i \in R$ and $j \in C$, it holds $\deg H_{i,j}(s) \leq 0$, which implies $\bar{c}_{i, j} \leq q_j - p_i$.
  Therefore $\prn{p, q}$ is feasible on $\D{\bar{A}}$.

  Next, we show that $\prn{p, q}$ is not optimal on $\D{\bar{A}}$.
  From~\eqref{eq:update_A_2}, the tight coefficient matrix $\bar{A}^\#$ of $\bar{A}(s)$ with respect to $\prn{p, q}$ is
  %
  \begin{align} \label{def:tilde_A_sharp}
    \bar{A}^\#
    =
    \begin{pmatrix} U & O \\ O & I \end{pmatrix} A^\#
    =
    \begin{pmatrix} \bar{Q}^\# \\ T^\# \end{pmatrix},
  \end{align}
  %
  where $\bar{Q}^\# = UQ^\#$.
  From Proposition~\ref{prop:trank_identity_min} and~\eqref{eq:rank_Qsharp}, it holds
  %
  \begin{align*}
    \trank \bar{A}^\#
    &= \min \set{\trank \bar{Q}^\#[R_Q, J] + \trank T^\#[R_T, J] + \card{C \setminus J}}[J \subseteq C] \\
    &\leq \trank \bar{Q}^\#[R_Q, J^*] + \trank T^\#[R_T, J^*] + \card{C \setminus J^*} \\
    &= \rank \bar{Q}^\#[R_Q, J^*] + \trank T^\#[R_T, J^*] + \card{C \setminus J^*}.
  \end{align*}
  %
  Now since $Q^\#[R_Q, J^*]$ and $\bar{Q}^\#[R_Q, J^*] = UQ^\#[R_Q, J^*]$ have the same rank, we obtain
  %
  \begin{align*}
    \trank \bar{A}^\#
    &\leq \rank Q^\#[R_Q, J^*] + \trank T^\#[R_T, J^*] + \card{C \setminus J^*} = \rank A^\#,
  \end{align*}
  %
  where the last equality comes from~\eqref{eq:Jstar}.
  In addition, since $\rank \bar{A}^\# = \rank A^\#$ from~\eqref{def:tilde_A_sharp}, we have $\trank \bar{A}^\# \leq \rank \bar{A}^\#$, which implies $\trank \bar{A}^\# = \rank \bar{A}^\# = \rank A^\#$.
  Furthermore, since $A(s)$ is not upper-tight, we have $\rank A^\# < n$ by Lemma~\ref{lem:tightness1}.
  Thus, $\trank A^\# = \rank A^\# < n$ holds.
  It then follows from Lemma~\ref{lem:tightness2} that $\prn{p, q}$ is not optimal on $\D{\bar{A}}$.
\end{proof}

From Lemma~\ref{lem:tilde_A_optimality} and the unimodularity of $U_Q(s)$, we obtain the following.

\begin{corollary}
  Let $A(s) = \binom{Q(s)}{T(s)}$ be an $n \times n$ nonsingular LM-polynomial matrix that is not upper-tight, and $\bar{A}(s)$ the LM-polynomial matrix defined in~\eqref{eq:update_A_2}.
  Then $\tdegdet{\bar{A}} \leq \tdegdet{A} - 1$ and $\degdet{A} = \degdet{\bar{A}}$ hold.
\end{corollary}

\subsection{Dual Updates}
\label{sec:dual_updates}

Let $\prn{p, q}$ be a feasible solution of $\D{\bar{A}}$.
We obtain an optimal solution of $\D{\bar{A}}$ by iterating the following procedure.

Let $\bar{A}^\#$ be the tight coefficient matrix of $\bar{A}(s)$ with respect to $\prn{p, q}$.
First we check $\trank \bar{A}^\# = n$ or not.
If it is, $\prn{p, q}$ is an optimal solution of $\D{\bar{A}}$ from Lemma~\ref{lem:tightness2} and we are done.
Otherwise, we construct a feasible solution $(p', q')$ of $\D{\bar{A}}$ such that the difference
\begin{align} \label{def:Delta}
  \Delta \defeq \sum_{j \in C} \prn{q'_j - q_j} - \sum_{i \in R} \prn{p'_i - p_i}
\end{align}
of the objective values is negative.
Let $G^\# = \prn{R, C; E^\#}$ be a bipartite graph defined by
\begin{align*}
  E^\#
  \defeq \set[\big]{\prn{i, j} \in R \times C}[\bar{A}^\#_{i,j} \neq 0]
  =      \set{\prn{i, j} \in \app{E}{\bar{A}}}[q_j - p_i = \bar{c}_{i,j}].
\end{align*}
Since $\prn{p, q}$ is not optimal, there is no perfect matching of $G^\#$.
Thus $G^\#$ has a vertex cover $S \subseteq R \cup C$ with $\card{S} < n$ by K\"{o}nig--Egerv\'{a}ry's theorem.
Using this $S$, we define $\prn{p', q'}$ as follows:
\begin{align} \label{eq:update_pq}
  p'_i = \begin{cases}
    p_i     & \prn{i \in R \cap S} \\
    p_i + 1 & \prn{i \in R \setminus S}
  \end{cases}
  , \quad
  q'_j = \begin{cases}
    q_j + 1 & \prn{j \in C \cap S} \\
    q_j     & \prn{j \in C \setminus S}
  \end{cases}
\end{align}
for $i \in R$ and $j \in C$.

\begin{lemma} \label{lem:dual_update}
  Let $\prn{p, q}$ be a feasible but not optimal solution of $\D{\bar{A}}$ and $\prn{p', q'}$ defined in~\eqref{eq:update_pq}.
  Then the difference $\Delta$ of the objective values in~\eqref{def:Delta} is negative, and $\prn{p', q'}$ is a feasible solution of $\D{\bar{A}}$.
\end{lemma}
\begin{proof}
  The difference of the objective values is $\Delta = \card{C \cap S} - \card{R \setminus S} = \card{S} - \card{R} < 0$.
  Next, we show the feasibility of $\prn{p', q'}$.
  For every $\prn{i, j} \in E(\bar{A})$, it holds $q_j - p_i \geq \bar{c}_{i,j}$ since $\prn{p, q}$ is feasible.
  If $i \in S$ or $j \in S$, it holds $\prn{q'_j - q_j} - \prn{p'_i - p_i} \geq 0$, which imply $q'_j - p'_i \geq q_j - p_i \geq \bar{c}_{i,j}$.
  If $i \notin W$ and $j \notin W$, then $\prn{i, j}$ is not an edge of $G^\#$ since $\prn{i, j}$ is not covered by $W$.
  Hence it holds $q_j - p_i > \bar{c}_{i,j}$, and thus $q'_j - p'_i = q_j - p_i - 1 \geq \bar{c}_{i,j}$.
\end{proof}

We update $\prn{p, q}$ to $\prn{p', q'}$, and go back to the optimality checking.
From Lemma~\ref{lem:dual_update}, it is guaranteed that $\prn{p, q}$ becomes an optimal solution of $\D{\bar{A}}$ by iterating the update process above.

\subsection{Complexity Analysis}
\label{sec:complexity_analysis}

This section is devoted to complexity analysis.
The dominating part in our algorithm is the matrix multiplications in~\eqref{eq:update_A_2}.

Let $A(s)$ be an $n \times n$ nonsingular LM-polynomial matrix and $A^\#$ the tight coefficient matrix with respect to an optimal solution $\prn{p, q}$ of $\D{A}$.
From the definition of $A^\#$, we can express $A(s)$ as
\begin{align} \label{eq:As_neg}
  A(s) =
    D_p^{-1}(s)
    \prn{A^\# + \sum_{k=1}^K s^{-k} V_k}
    D_q(s)
\end{align}
for some $K$ matrices $V_1, V_2, \ldots, V_K$ with $V_K \neq O$.
By~\eqref{eq:update_A_2} and~\eqref{eq:As_neg}, we have
\begin{align*}
  \bar{A}(s) =
    D^{-1}_p(s)
    \begin{pmatrix} U & O \\ O & I \end{pmatrix}
    \prn{A^\# + \sum_{k=1}^K s^{-k} V_k}
    D_q(s).
\end{align*}
Therefore, we can compute $\bar{A}(s)$ by performing $K+1$ constant matrix multiplications.

By $V_K \neq O$, there exist $i \in R$ and $j \in C$ such that the $\prn{i, j}$ entry in $V_K$ is nonzero.
Then the degree of the corresponding term in $A_{i,j}(s)$ is equal to $q_j - p_i - K$.
Since $A_{i,j}(s)$ is a polynomial, we have $q_j - p_i - K \geq 0$, which implies $K \leq q_j - p_i \leq q_j$.
The following lemma bounds $p_i$ and $q_j$ at any iteration of our algorithm.

\begin{lemma} \label{lem:bound_l}
  During the algorithm, the values $p_i$ and $q_j$ are at most $2 ln$ for $i \in R$ and $j \in C$, where $l$ is the maximum degree of an entry in $A(s)$.
\end{lemma}
\begin{proof}
  From Lemma~\ref{lem:bound_pq}, the initial values of $p_i$ and $q_j$ are bounded by $ln$.
  In every update of $\prn{p, q}$, the values $p_i$ and $q_j$ increases by at most one from the update rule~\eqref{eq:update_pq}.
  In addition, $\prn{p, q}$ is updated at most $\tdegdet{A} - \degdet{A} \leq ln$ times because the objective value $\sum_{j \in C} q_j - \sum_{i \in R} p_i$ of the dual problem decreases by at least one in every update.
  Therefore, at any iteration of the algorithm, it holds $p_i, q_j \leq ln + \tdegdet {A} \leq ln + ln = 2ln$.
\end{proof}

The time complexity of our algorithm is as follows.
\begin{theorem} \label{thm:complexity_tightness}
  Let $A(s)$ be an $n \times n$ nonsingular LM-polynomial matrix and $l$ the maximum degree of an entry in $A(s)$.
  Then Algorithm for Tightness runs in $\Order{l^2 n^{\omega + 2}}$ time, where $2 < \omega \leq 3$ is the matrix multiplication exponent.
\end{theorem}

\begin{proof}
  Phase~1 can be done in $\Order{n^3}$ time by the Hungarian method~\cite{Kuhn1955} and shortest path algorithms such as the Bellman--Ford algorithm.
  Consider the time complexity in every iteration of Phases~2 and~3.
  In Phase~2, the nonsingularity of the tight coefficient matrix $A^\#$ can be checked via the rank identity~\eqref{eq:rank_identity_min}.
  Thus an efficient way is to obtain a minimizer $J^*$ of~\eqref{eq:rank_identity_min} before Phase 2, and then check the nonsingularity of $A^\#$ by~\eqref{eq:rank_identity_min}.
  The minimizer $J^*$ can be found from a residual graph constructed by an augmenting path type algorithm~\cite{Murota1987}, which runs in $\Order{n^3 \log n}$ time~\cite{Cunningham1986}.
  The computation of $\bar{A}(s)$ in Phase~3 can be done in $\Order{N n^\omega} = \mathrm{O}\prn[\big]{\max_{j \in C} q_j n^\omega} = \Order{l n^{\omega + 1}}$ time from Lemma~\ref{lem:bound_l}, where $\prn{p, q}$ is a dual optimal solution of $\D{A}$ and $N$ is in~\eqref{eq:As_neg}.
  In addition, since the number of iterations of Phases~2 and~3 is at most $\tdegdet{A} - \degdet{A} \leq ln$, the running time in Phases~2 and~3 is $\Order{l^2 n^{\omega+2}}$.
  Finally, the updates of $\prn{p, q}$ run in $\Order{ln^4}$ time: $\prn{p, q}$ is updated at most $\tdegdet{A} \leq ln$ times, and in every update, we can find a vertex cover in $\Order{n^3}$ time by Ford--Fulkerson's algorithm.
  Thus the total running time is $\Order{l^2 n^{\omega+2}}$.
\end{proof}

\begin{theorem} \label{thm:complexity_all}
  For a DAE~\eqref{eq:laplaced_mixed_dae} with $n \times n$ nonsingular mixed polynomial matrix $A(s)$, our algorithm returns an equivalent DAE of index zero or one in $\Order{l^2 n^{\omega+2}}$ time, where $2 < \omega \leq 3$ is the matrix multiplication exponent and $l$ is the maximum degree of entries in $A(s)$.
\end{theorem}

\begin{proof}
  First we convert the DAE into an equivalent DAE with LM-polynomial matrix $\LM{A}(s)$ of size $2n \times 2n$.
  Note that the maximum degree of an entry in $\LM{A}(s)$ is equal to $l$ by~\eqref{eq:LM_dae}.
  Hence it follows from Theorem~\ref{thm:complexity_tightness} that Algorithm for Tightness for $\LM{A}(s)$ runs in $\Order{l^2 n^{\omega+2}}$ time.
  The resulting DAE has a coefficient matrix such that the maximum degree of an entry is at most $4ln$, because it holds that
  %
  \begin{align*}
    \deg \LM{A}_{i,j}(s) \leq q_j - p_i \leq q_j \leq 4ln
  \end{align*}
  %
  with a feasible solution $\prn{p, q}$ of $\D{\LM{A}}$, where the last inequality is due to Lemma~\ref{lem:bound_l}.

  Next we analyze the complexity of the MS-algorithm described in Section~\ref{sec:ms_alg}.
  In Step~1, we can reuse a dual optimal solution $\prn{p, q}$ obtained at the termination of Algorithm for Tightness, or compute a new $\prn{\tilde{p}, \tilde{q}}$ such that $\tilde{p}_i \leq p_i$ for $i \in R$ to decrease the number of dummy variables, in $\Order{n^3}$ time.
  The nonsingularity of the corresponding tight coefficient matrix can be verified by solving an independent matching problem in $\Order{n^3 \log n}$ time~\cite{Cunningham1986, Murota1987}.
  Step~2 runs in $\Order{n^4 \log n}$ time since we solve independent matching problems at most $2n$ times.
  We now consider the resultant DAE returned in Step~4.
  The number of original (non-dummy) variables is $2n$, and from Lemma~\ref{lem:bound_l}, the orders of their derivatives are at most
  %
  \begin{align*}
    4ln + \max_{i \in R} p_i \leq 4ln + 4ln = \Order{ln}.
  \end{align*}
  %
  In contrast, the number of dummy variables is $\sum_{i \in R} p_i = \Order{ln^2}$, and there is no derivative of dummy variables in the resultant DAE.
  Therefore, the number of terms in the resultant DAE is $\Order{l n^2}$, and thus Step~4 runs in $\Order{ln^2}$ time.
  Therefore the MS-algorithm costs $\Order{n^4 \log n + l n^2}$ time.

  Since the bottleneck in the entire algorithm is Algorithm for Tightness, the total running time of our algorithm is $\Order{l^2 n^{\omega+2}}$.
\end{proof}

\section{Exploiting Dimensional Consistency}
\label{sec:exploiting_dimensional_consistency}

\subsection{Dimensional Consistency}

The principle of dimensional homogeneity claims that any equation describing a physical phenomenon must be consistent with respect to physical dimensions.
In order to reflect the dimensional consistency in conservation laws of dynamical systems, Murota~\cite{Murota1985b} introduced a class of mixed polynomial matrices $A(s) = Q(s) + T(s)$ that satisfy the following condition.
\begin{enumerate}[label={(MP-DC)}]
  \item[(MP-DC)] Every nonvanishing subdeterminant of $Q(s)$ is a monomial in $s$.
\end{enumerate}
A mixed polynomial matrix satisfying (MP-DC) is said to be \emph{dimensionally consistent}.
We abbreviate a dimensionally consistent mixed polynomial matrix and a dimensionally consistent LM-polynomial matrix to a \emph{DCM-polynomial matrix} and a \emph{DCLM-polynomial matrix}, respectively.
It is known that an $m \times n$ polynomial matrix $Q(s)$ satisfies (MP-DC) if and only if $Q(s)$ is written as
\begin{align} \label{eq:dimensional_consistency_1}
  Q(s) =
  \diag{s^{-\lambda_1}, \ldots, s^{-\lambda_{m}}}
  Q(1)
  \diag{s^{\mu_1}, \ldots, s^{\mu_n}}
\end{align}
for some integers $\lambda_1, \ldots, \lambda_{m}$ and $\mu_1, \ldots, \mu_n$.

\subsection{Improved Algorithm}
\label{sec:improved_algorithm}

This section improves the matrix modification procedure in Phase~3 for DCLM-polynomial matrices preserving their dimensional consistency.

Let $A(s) = \binom{Q(s)}{T(s)}$ be a DCLM-polynomial matrix with $R_Q = \row{Q}$, $R_T = \row{T}$ and $C = \col{A}$.
Let $A^\# = \binom{Q^\#}{T^\#}$ denote the tight coefficient matrix with respect to an optimal solution $\prn{p, q}$ of $\D{A}$.
For an integer $k \in \setZ$, let
\begin{align} \label{def:Rk_Ck}
  R_k = \set{i \in R_Q}[p_i - \lambda_i = k],
  \quad
  C_k = \set{j \in C}[q_j - \mu_j = k].
\end{align}
If $Q_{i,j}^\# \neq 0$, then we have $c_{i,j} = q_j - p_i$ from the definition of $A^\#$ and $c_{i,j} = \mu_j - \lambda_i$ by~\eqref{eq:dimensional_consistency_1}.
Hence $q_j - p_i = \mu_j - \lambda_i$ follows, which implies $i \in R_k$ if and only if $j \in C_k$.
Thus, it holds $A^\#[R_h, C_k]=O$ for distinct $h,k \in \setZ$.
Namely, $Q^\#$ forms a block diagonal matrix as
\begin{align*}
  Q^\# = \begin{blockarray}{ccccccc}
    & \cdots & C_{-1} & C_0 & C_1 & C_2 & \cdots \\
    \begin{block}{c(cccccc)}
      \vdots & \ddots \\
      R_{-1} && Q^\#_{-1} \\
      R_{0}  &&& Q^\#_0 \\
      R_{1}  &&&& Q^\#_1 \\
      R_{2}  &&&&& Q^\#_2 \\
      \vdots &&&&&& \ddots \\
    \end{block}
  \end{blockarray}
  \,\, ,
\end{align*}
where $Q^\#_k = Q^\#[R_k, C_k]$ for $k \in \setZ$, and empty blocks indicate zero submatrices.

Let $J^* \subseteq C$ be a minimizer of the rank identity~\eqref{eq:rank_identity_min} for $A^\#$.
Sorting rows in ascending order of $p$, the matrix modification process described in Section~\ref{sec:matrix_modification} finds a nonsingular upper-triangular matrix $U$ such that
\begin{align} \label{eq:UQsharp}
  \rank UQ^\#[R_Q, J^*] = \trank UQ^\#[R_Q, J^*].
\end{align}
For a DCLM-polynomial matrix, supposing that rows in $R_k$ is sorted in ascending order of $p$, we find a nonsingular upper-triangular matrix $U_k$ such that
\begin{align*}
  \rank U_k Q_k^\#[R_k, C_k \cap J^*] = \trank U_k Q_k^\#[R_k, C_k \cap J^*]
\end{align*}
for $k \in \setZ$.
Then $U = \blockdiag{\ldots, U_{-1}, U_0, U_1, U_2, \ldots}$ satisfies~\eqref{eq:UQsharp}, where $\blockdiag{B_1, B_2, \ldots, B_N}$ is a block diagonal matrix of diagonal blocks $B_1, B_2, \ldots, B_N$.

For $k \in \setZ$, let $P_k(s)$ be a diagonal polynomial matrix with $\row{P_k} = \col{P_k} = R_k$ whose $(i, i)$ entry is $s^{p_i}$ for each $i \in R_k$.
Let $D_p(s) = \blockdiag{\ldots, P_{-1}(s), P_0(s), P_1(s), P_2(s), \ldots}$.
Now the unimodular matrix $U_Q(s)$ defined in~\eqref{def:U_tilde} can be written as
\begin{align}
  U_Q(s)
  &= D_p^{-1}(s) \blockdiag{\ldots, U_{-1}, U_0, U_1, U_2, \ldots} D_p(s) \nonumber \\
  &= D_p^{-1}(s) \blockdiag{\ldots, U_{-1} P_{-1}(s), U_0 P_0(s),  U_1 P_1(s), U_2 P_2(s), \ldots}. \label{eq:DC_UQ}
\end{align}
Then we update $A(s)$ into $\bar{A}(s) = \binom{U_Q(s)Q(s)}{T(s)}$ as written in~\eqref{eq:update_A_2}.

\begin{lemma}
  Let $A(s) = \binom{Q(s)}{T(s)}$ be an $n \times n$ DCLM-polynomial matrix.
  Then $\bar{A}(s) = \binom{U_Q(s)Q(s)}{T(s)}$ is also dimensionally consistent.
\end{lemma}
\begin{proof}
  Let $\lambda_1, \ldots, \lambda_{m_Q}$ and $\mu_1, \ldots, \mu_n$ defined in~\eqref{eq:dimensional_consistency_1} for $A(s)$, where $m_Q = \card{\row{Q}}$.
  For $k \in \setZ$, let $R_k$ and $C_k$ defined in~\eqref{def:Rk_Ck}, and let $\Lambda_k(s)$ denote a diagonal polynomial matrix with $\row{\Lambda_k} = \col{\Lambda_k} = R_k$ whose $(i, i)$ entry is $s^{\lambda_i}$ for each $i \in R_k$, and $D_\mu(s) = \diag{s^{\mu_1}, \ldots, s^{\mu_n}}$.
  Then the condition~\eqref{eq:dimensional_consistency_1} for dimensional consistency is written as
  %
  \begin{align} \label{eq:DC_Q}
    Q(s) = \blockdiag{\ldots, \Lambda_{-1}^{-1}(s), \Lambda_0^{-1}(s), \Lambda_1^{-1}(s), \Lambda_2^{-1}(s), \ldots} Q(1) D_\mu(s).
  \end{align}
  %
  Combining~\eqref{eq:DC_UQ} and~\eqref{eq:DC_Q}, we obtain
  %
  \begin{align}
    U_Q(s) Q(s)
    &= P^{-1}(s) \blockdiag{\ldots, U_{-1} P_{-1}(s) \Lambda_{-1}^{-1}(s), U_0 P_0(s) \Lambda_0^{-1}(s),  U_1 P_1(s) \Lambda_1^{-1}(s), \ldots} Q(1) D_\mu(s) \nonumber \\
    &= P^{-1}(s) \blockdiag{\ldots, s^{-1} U_{-1} , U_0,  s U_1, s^2 U_2, \ldots} Q(1) D_\mu(s) \nonumber \\
    &= \blockdiag{\ldots, s^{-1} P_{-1}^{-1}(s), P_0^{-1}(s), sP_1^{-1}(s), s^2P_2^{-1}(s), \ldots} UQ(1)D_\mu(s) \label{eq:DC_tilde_Q},
  \end{align}
  %
  where we used $P_k(s) \Lambda_k^{-1}(s) = s^k I$ for $k \in \setZ$.
  From~\eqref{eq:DC_tilde_Q}, $\bar{A}(s)$ is also dimensionally consistent.
\end{proof}

\subsection{Complexity Analysis}
For a DCLM-polynomial matrix $A(s)$, we can compute $\bar{A}(s) = U(s) A(s)$ only by one constant matrix multiplication $U Q(1)$ from~\eqref{eq:DC_tilde_Q}, whereas a general LM-polynomial matrix needs $\Order{ln}$ multiplications.
This improves the total running time as follows.

\begin{theorem} \label{thm:DC_complexity_tightness}
  Let $A(s)$ be an $n \times n$ nonsingular DCLM-polynomial matrix and $l$ the maximum degree of an entry in $A(s)$.
  Then Algorithm for Tightness runs in $\Order{l n^4 \log n}$ time.
\end{theorem}
\begin{proof}
  For each iteration of Phases~2 and~3, the computation of $\bar{A}(s)$ in Phase~3 can be done in $\Order{n^\omega}$ time, where $2 < \omega \leq 3$ is the matrix multiplication exponent.
  The most expensive part is the nonsingularity checking for a tight coefficient matrix in Phase~2, which requires $\Order{n^3 \log n}$ time~\cite{Cunningham1986, Murota1987}.
  Since the number of iterations of Phases~2 and~3 is at most $\tdegdet{A} - \degdet{A} \leq ln$, the running time of Phases~2 and~3 is $\Order{l n^4 \log n}$.
  We can check that other computations run in $\Order{ln^4 \log n}$ time as in the proof of Theorem~\ref{thm:complexity_tightness}.
\end{proof}

\begin{theorem} \label{thm:DC_complexity_all}
  For a DAE~\eqref{eq:laplaced_mixed_dae} with $n \times n$ nonsingular DCM-polynomial coefficient matrix $A(s)$, our algorithm returns an equivalent DAE of index zero or one in $\Order{l n^4 \log n}$ time, where $l$ is the maximum degree of entries in $A(s)$.
\end{theorem}
\begin{proof}
  We can easily check that the coefficient LM-polynomial matrix of the augmented DAE described in Section~\ref{sec:mixed_to_LM} is also dimensionally consistent.
  Algorithm for Tightness runs in $\Order{ln^4 \log n}$ time from Theorem~\ref{thm:DC_complexity_tightness}.
  In addition, the MS-algorithm runs in $\Order{n^4 \log n + ln^2}$ time as discussed in the proof of Theorem~\ref{thm:complexity_all}.
  Thus the total running time is $\Order{ln^4 \log n}$.
\end{proof}

\section{Examples}
\label{sec:example}

We give two examples below.
The first example is a simple index-4 DAE and the second example is a DAE representing an electrical network.
Throughout the execution of our algorithm, it is emphasized that: (i) we only use combinatorial operations and numerical calculations over rational numbers (especially over integers in the following examples), and (ii) we do not reference nominal values of physical quantities.

\subsection{Example of High-index DAE}
The first example is the following index-4 DAE
\begin{align} \label{eq:example1}
  \left\{\begin{aligned}
    \ddot{x}_1 - \dot{x}_1 + \ddot{x}_2 - \dot{x}_2 + x_4 &= f_1(t), \\
    \ddot{x}_1 + \ddot{x}_2 + x_3 &= f_2(t), \\
    \alpha_1 x_2 + \alpha_2 \ddot{x}_3 + \alpha_3 \dot{x}_4 &= f_3(t), \\
    \alpha_4 x_3 + \alpha_5 \dot{x}_4 &= f_4(t),
  \end{aligned}\right.
\end{align}
with independent parameters $\alpha_1, \ldots, \alpha_5$ and smooth functions $f_1, \ldots, f_4$.
The coefficient matrix $A(s) = \binom{Q(s)}{T(s)}$ corresponding to~\eqref{eq:example1} is an LM-polynomial matrix given by
\begin{align} \label{eq:example1_A}
  A(s) = \begin{pmatrix}
    s^2-s & s^2-s    &              & 1 \\
    s^2   & s^2      & 1            &   \\
          & \alpha_1 & \alpha_2 s^2 & \alpha_3 s \\
          &          & \alpha_4     & \alpha_5 s
  \end{pmatrix},
\end{align}
where empty cells indicate zero.
The row sets $R_Q$ of $Q(s)$ and $R_T$ of $T(s)$ correspond to the first and last two rows in $A(s)$, respectively.
Since $\degdet{A} = \deg \prn{-\alpha_1 \alpha_5 s^3 -\alpha_1 \alpha_4 s^2 + \alpha_1 \alpha_5 s^2} = 3$ and $\tdegdet{A} = 7$, the MS-algorithm is not applicable to the DAE, which is shown in our algorithm.

Let us apply our algorithm to~\eqref{eq:example1_A}.
First, we find a dual optimal solution $p = \prn{0, 0, 0, 0}$ and $q = \prn{2, 2, 2, 1}$.
The corresponding tight coefficient matrix $A^\# = \binom{Q^\#}{T^\#}$ is
\begin{align*}
  A^\# = \begin{pmatrix}
    1 & 1 &   & \\
    1 & 1 &   & \\
      &   & \alpha_2 & \alpha_3 \\
      &   &          & \alpha_5 \\
  \end{pmatrix}.
\end{align*}
A minimizer $J^*$ of~\eqref{eq:rank_identity_min} for $A^\#$ is a set of the first two columns as follows:
\newcommand\overmat[3]{%
  \BAmulticolumn{#1}{c}{%
    \overbrace{%
      \hphantom{%
        \begin{matrix}#2\end{matrix}%
      }%
    }^{#3}%
  }%
}
\newcommand\rightmat[2]{\rdelim\}{#1}{0em}[$\scriptstyle{#2}$]}
\begin{align*}
  A^\# = \begin{blockarray}{ccccc}
    \overmat{2}{1 & 1}{J^*} & \overmat{2}{\alpha_2 & \alpha_3}{C \setminus J^*}
    \\
    \begin{block}{(cccc)c}
      1 & 1 &          &          & \rightmat{2}{R_Q} \\
      1 & 1 &          &          & \\
        &   & \alpha_2 & \alpha_3 & \rightmat{2}{R_T} \\
        &   &          & \alpha_5 & \\
    \end{block}
  \end{blockarray}
  \quad \quad.
\end{align*}
Then we can check that $\rank A^\# = Q^\#[R_Q, J^*] + T^\#[R_T, J^*] + \card{C \setminus J^*} = 1 + 0 + 2 = 3 < 4$, which implies that $A(s)$ is not upper-tight.
We convert $Q^\#[R_Q, J^*] = \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}$ by the backward elimination into
\begin{align*}
  \bar{Q}^\#[R_Q, J^*] = U Q^\#[R_Q, J^*] = \begin{pmatrix} & \\ 1 & 1 \end{pmatrix},
\end{align*}
where $U = \begin{pmatrix} 1 & -1 \\ & 1 \end{pmatrix}$.
Using $U_Q(s) = U$, the LM-polynomial matrix $A(s)$ is modified to
\begin{align*}
  A'(s) =
  \begin{pmatrix}
    1 & -1 &   &  \\
      &  1 &   &  \\
      &    & 1 &  \\
      &    &   & 1
  \end{pmatrix}
  A(s)
  =
  \begin{pmatrix}
    -s  & -s       & -1           & 1 \\
    s^2 & s^2      &  1           &   \\
        & \alpha_1 & \alpha_2 s^2 & \alpha_3 s \\
        &          & \alpha_4     & \alpha_5 s
  \end{pmatrix}.
\end{align*}
The dual solution is updated to $p' = \prn{1, 0, 0, 1}$ and $q' = \prn{2, 2, 2, 2}$, and the corresponding tight coefficient matrix $A'^\# = \binom{Q'^\#}{T'^\#}$ of $A'(s)$ is
\begin{align*}
  A'^\# = \begin{pmatrix}
    -1 & -1 &    & \\
     1 &  1 &    & \\
       &    & \alpha_2 & \\
       &    &          & \alpha_5 \\
  \end{pmatrix}.
\end{align*}
The minimizer $J^*$ that we used above also minimizes the right-hand side of the rank identity~\eqref{eq:rank_identity_min} for $A'^\#$.
Since $A'^\#$ is still singular, we go on the modification.
Noting the order of rows, we transform $Q'^\#[R_Q, J^*] = \begin{pmatrix} -1 & -1 \\ 1 & 1 \end{pmatrix}$ by $U' = \begin{pmatrix} 1 & \\ 1 & 1 \end{pmatrix}$ into
\begin{align*}
  \bar{Q}'^\#[R_Q, J^*] = U' Q'^\#[R_Q, J^*] = \begin{pmatrix} -1 & -1 \\ & \end{pmatrix}.
\end{align*}
We have $U_Q'(s) = \diag{s^{-1}, 1} U' \diag{s, 1} = \begin{pmatrix} 1 & \\ s & 1 \end{pmatrix}$, and modify $A'(s)$ to
\begin{align*}
  A''(s) =
  \begin{pmatrix}
    1 &   &   &  \\
    s & 1 &   &  \\
      &   & 1 &  \\
      &   &   & 1
  \end{pmatrix}
  A'(s)
  =
  \begin{pmatrix}
    -s & -s       & -1           & 1   \\
       &          & -s+1         & s   \\
       & \alpha_1 & \alpha_2 s^2 & \alpha_3 s \\
       &          & \alpha_4     & \alpha_5 s
  \end{pmatrix}.
\end{align*}
The dual solution is updated to $p'' = \prn{1, 3, 2, 3}$ and $q'' = \prn{2, 2, 4, 4}$.
Our algorithm halts at this point since $A''(s)$ is upper-tight, which can be checked through the nonsingularity of the tight coefficient matrix $A''^\#$ again.
Now $\degdet{A}$ is computed as $\degdet{A} = \degdet{A''} = \tdegdet{A''} = 3$.
The resulting DAE is
\begin{align} \label{eq:example1_result}
  \left\{\begin{aligned}
    -\dot{x}_1 -\dot{x}_2 -x_3 + x_4 &= f_1(t) - f_2(t), \\
    -\dot{x}_3 + x_3 + \dot{x}_4 &= \dot{f}_1(t) - \dot{f}_2(t) + f_2(t),\\
     \alpha_1 x_2 + \alpha_2 \ddot{x}_3 + \alpha_3 \dot{x}_4 &= f_3(t), \\
     \alpha_4 x_3 + \alpha_5 \dot{x}_4 &= f_4(t),
  \end{aligned}\right.
\end{align}
which is index two.

An index-1 DAE is obtained by applying the MS-algorithm to the DAE~\eqref{eq:example1_result}.
Instead of $\prn{p'', q''}$, we now use an optimal solution $\tilde{p} = \prn{0,2,1,2}$ and $\tilde{q} = \prn{1,1,3,3}$ of $\D{A''}$ in the MS-algorithm in order to decrease the number of dummy variables and equations.
Then the MS-algorithm outputs an index-1 DAE
\begin{align*}
  \left\{\begin{aligned}
    -\dot{x}_1 -z^{[1]}_2 -x_3 + x_4 &= f_1(t) - f_2(t), \\
    -\dot{x}_3 + x_3 + \dot{x}_4 &= \dot{f}_1(t) - \dot{f}_2(t) + f_2(t),\\
    -z^{[2]}_3 + x_3 + z^{[2]}_4 &= \ddot{f}_1(t) - \ddot{f}_2(t) + \dot{f}_2(t),\\
    -z^{[3]}_3 + x_3 + z^{[3]}_4 &= \dddot{f}_1(t) - \dddot{f}_2(t) + \ddot{f}_2(t),\\
     \alpha_1 x_2 + \alpha_2 x_3 + \alpha_3 \dot{x}_4 &= f_3(t), \\
     \alpha_1 z^{[1]}_2 + \alpha_2 \dot{x}_3 + \alpha_3 z^{[2]}_4 &= \dot{f}_3(t), \\
     \alpha_4 x_3 + \alpha_5 \dot{x}_4 &= f_4(t), \\
     \alpha_4 \dot{x}_3 + \alpha_5 z^{[2]}_4 &= \dot{f}_4(t), \\
     \alpha_4 z^{[2]}_3 + \alpha_5 z^{[3]}_4 &= \ddot{f}_4(t), \\
  \end{aligned}\right.
\end{align*}
where $z^{[1]}_2$, $z^{[2]}_3$, $z^{[3]}_3$, $z^{[2]}_4$ and $z^{[3]}_4$ are dummy variables corresponding to $\dot{x}_2$, $\ddot{x}_3$, $\dddot{x}_3$, $\ddot{x}_4$ and $\dddot{x}_4$, respectively.

\subsection{Example of Electrical Network}
\label{sec:example_2}

The next example is a DAE representing an electrical network illustrated in Figure~\ref{fig:RLC}, given in~\cite[Section 1.1]{Murota2000}.
The network consists of a voltage source of time-varying voltage $V(t)$, two resistances $R_1$ and $R_2$, an inductor $L$ and a capacitor $C$.
State variables of this network is currents $i_1, \ldots, i_5$ shown in Figure~\ref{fig:RLC}, and voltages $v_1, \ldots, v_5$ across branches carrying the currents $i_1, \ldots, i_5$, respectively.

\begin{figure}[htbp]
  \centering
  \begin{circuitikz}[]
    \draw (0,-2)
      to[short]  (0,-1)
      to[sV=$V(t)$] (0,1)
      to[short, i=$i_5$]  (0,2)
      to[short]  (2,2)
      to[resistor=$R_1$, i=$i_1$] (2,0)
      to[resistor=$R_2$, i=$i_2$] (4,0);
    \draw (2,2)
      to[short] (4,2)
      to[C=$C$, i=$i_4$] (4,0)
      to[short] (4,-2)
      to[short] (2,-2);
    \draw (2,0)
      to[L=$L$, i=$i_3$] (2,-2)
      to[short] (0,-2);
  \end{circuitikz}
  \caption{%
    Simple RLC network.
  }
  \label{fig:RLC}
\end{figure}
The Laplace transform of an index-2 DAE representing this network is given by
\begin{align} \label{eq:example2_eq}
  \prn{\begin{array}{rrrrr|rrrrr}
    -1  &     &     & -1 &  1 &    &    &    &    &    \\
        &   1 &   1 &  1 & -1 &    &    &    &    &    \\
    \hline
        &     &     &    &    &  1 &    &  1 &    & -1 \\
        &     &     &    &    & -1 & -1 &    &  1 &    \\
        &     &     &    &    &    &  1 & -1 &    &    \\
    \hline
    R_1 &     &     &    &    & -1 &    &    &    &    \\
        & R_2 &     &    &    &    & -1 &    &    &    \\
        &     &  sL &    &    &    &    & -1 &    &    \\
        &     &     & -1 &    &    &    &    & sC &    \\
        &     &     &    &    &    &    &    &    &  1
  \end{array}}
  \begin{pmatrix}
    \tilde{i}_1 \\
    \tilde{i}_2 \\
    \tilde{i}_3 \\
    \tilde{i}_4 \\
    \tilde{i}_5 \\
    \tilde{v}_1 \\
    \tilde{v}_2 \\
    \tilde{v}_3 \\
    \tilde{v}_4 \\
    \tilde{v}_5
  \end{pmatrix}
  =
  \begin{pmatrix}
    0 \\
    0 \\
    0 \\
    0 \\
    0 \\
    0 \\
    0 \\
    0 \\
    0 \\
    \tilde{V}(s)
  \end{pmatrix},
\end{align}
where $\tilde{i}_1, \ldots, \tilde{i}_5$, $\tilde{v}_1, \ldots, \tilde{v}_5$ and $\tilde{V}(s)$ are the Laplace transforms of $i_1, \ldots, i_5$, $v_1, \ldots, v_5$ and $V(t)$, respectively (we assumed that all state variables and their derivatives were equal to zero at $t = 0$ for simplicity).
In this system~\eqref{eq:example2_eq}, the first two equations come from Kirchhoff's current law (KCL), and the following three equations come from Kirchhoff's voltage law (KVL).
The last five equations represent the element characteristics (constitutive equations).
The coefficient matrix in~\eqref{eq:example2_eq} is naturally regarded as a mixed polynomial matrix with independent parameters $R_1$, $R_2$, $L$ and $C$ since nominal values of the parameters are supposed to be inaccurate.

Let $A(s)$ denote the coefficient matrix in~\eqref{eq:example2_eq}.
Since $A(s)$ is not LM-polynomial, it seems that we cannot directly apply our algorithm to $A(s)$.
However, since each of the last five rows in $A(s)$ do not contain two or more accurate constants, we can convert $A(s)$ into an LM-polynomial matrix by multiplying an independent parameter to each of the rows.
In addition, by the same logic to Lemma~\ref{lem:B}, our algorithm works without actually multiplying the independent parameters by regarding nonzero entries in the last five rows as independent parameters.
Thus we see $A(s)$ as an LM-polynomial matrix $A(s) = \binom{Q(s)}{T(s)}$, where $Q(s)$ and $T(s)$ correspond to the first and last five rows in $A(s)$, respectively.
In addition, $A(s)$ meets the condition~\eqref{eq:dimensional_consistency_1} for DCLM-polynomial matrices with $\lambda = \prn{0, 0, 0, 0, 0}$ and $\mu = \prn{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}$.

We are now ready for applying our algorithm to $A(s)$.
In Phase 1, a dual optimal solution is obtained as $p = \prn{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}$ and $q = \prn{0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0}$, which implies that $\tdegdet{A} = 2$.
The corresponding tight coefficient matrix $A^\# = \binom{Q^\#}{T^\#}$ is given by
\begin{align*}
  A^\# =
  \prn{\begin{array}{rrrrr|rrrrr}
    -1  &     &     & -1 &  1 &    &    &    &    &    \\
        &   1 &     &  1 & -1 &    &    &    &    &    \\
    \hline
        &     &     &    &    &  1 &    &  1 &    & -1 \\
        &     &     &    &    & -1 & -1 &    &    &    \\
        &     &     &    &    &    &  1 & -1 &    &    \\
    \hline
    R_1 &     &     &    &    & -1 &    &    &    &    \\
        & R_2 &     &    &    &    & -1 &    &    &    \\
        &     &   L &    &    &    &    & -1 &    &    \\
        &     &     & -1 &    &    &    &    &  C &    \\
        &     &     &    &    &    &    &    &    &  1
  \end{array}}
  \begin{array}{@{}l}
    \rdelim\}{5}{0em}[$\scriptstyle{R_Q}$]  \\
    \\
    \\
    \\
    \\
    \rdelim\}{5}{0em}[$\scriptstyle{R_T}$]  \\
    \\
    \\
    \\
    \\
  \end{array}
  \quad .
\end{align*}
A minimizer $J^*$ of the rank identity~\eqref{eq:rank_identity_min} for $A^\#$ is the set of nine columns other than the rightmost column corresponding to the variable $\tilde{v}_5$.
Thus we can check
\begin{align*}
  \rank A^\# = Q^\#[R_Q, J^*] + T^\#[R_T, J^*] + \card{C \setminus J^*} = 4 + 4 + 1 = 9 < 10,
\end{align*}
which imply that $A(s)$ is not upper-tight.
We proceed to the matrix modification process for DCLM-polynomial matrices that we described in Section~\ref{sec:improved_algorithm}.

The row set $R_k$ and the column set $C_k$ for $k \in \setZ$ defined in~\eqref{def:Rk_Ck} is the following:
\begin{align*}
  \begin{blockarray}{ccc@{}c@{}ccccc@{}c@{\,}cc}
    &\overmat{9}{-1&-1&-1&-1&-1&-1-1&-1&-1}{J^*} \\
    &\overmat{2}{-1&-1}{C_0} & \overbrace{}^{C_1} & \overmat{5}{-1&-1&-1&-1&-1}{C_0} & \overbrace{}^{C_1} & \overbrace{}^{C_0}\\
    \begin{block}{c(cc@{}c@{}ccccc@{}c@{\,}c)c}
      &-1  &     &     & -1 &  1 &    &    &    &   &    & \rightmat{5}{R_0 .} \\
      &    &   1 &     &  1 & -1 &    &    &    &   &    \\
      Q^\# = &    &     &     &    &    &  1 &    &  1 &   & -1 \\
      &    &     &     &    &    & -1 & -1 &    &   &    \\
      &    &     &     &    &    &    &  1 & -1 &   &    \\
    \end{block}
  \end{blockarray}
\end{align*}
Now $Q^\#$ can be seen as a block diagonal matrix consisting of one diagonal block $Q^\#_0 = Q^\#[R_0, C_0]$ by $Q^\#[R_0, C_1] = O$.
We transform
\begin{align*}
  Q_0^\#[R_0, C_0 \cap J^*]
  =
  \prn{\begin{array}{rrrrrrr}
    -1  &     &  -1 &  1 &    &    &     \\
        &   1 &   1 & -1 &    &    &     \\
        &     &     &    &  1 &    &   1 \\
        &     &     &    & -1 & -1 &     \\
        &     &     &    &    &  1 &  -1 \\
  \end{array}}
\end{align*}
into
\begin{align*}
  U Q_0^\#[R_0, C_0 \cap J^*]
  =
  \prn{\begin{array}{rrrrrrr}
    -1  &     &  -1 &  1 &    &    &     \\
        &   1 &   1 & -1 &    &    &     \\
        &     &     &    &    &    &     \\
        &     &     &    & -1 & -1 &     \\
        &     &     &    &    &  1 &  -1 \\
  \end{array}},
\end{align*}
where $U = \prn{\begin{array}{ccccc} 1 &&&& \\ & 1 &&& \\ && 1 & 1 & 1 \\ &&& 1 & \\ &&&& 1 \end{array}}$.
Using $U_Q(s) = U$, we modify $A(s)$ to
\begin{align*}
  A'(s) =
  \prn{\begin{array}{rrrrr|rrrrr}
    -1  &     &     & -1 &  1 &    &    &    &    &    \\
        &   1 &   1 &  1 & -1 &    &    &    &    &    \\
    \hline
        &     &     &    &    &    &    &    &  1 & -1 \\
        &     &     &    &    & -1 & -1 &    &  1 &    \\
        &     &     &    &    &    &  1 & -1 &    &    \\
    \hline
    R_1 &     &     &    &    & -1 &    &    &    &    \\
        & R_2 &     &    &    &    & -1 &    &    &    \\
        &     &  Ls &    &    &    &    & -1 &    &    \\
        &     &     & -1 &    &    &    &    & Cs &    \\
        &     &     &    &    &    &    &    &    &  1
  \end{array}},
\end{align*}
where the third row is different between $A(s)$ and $A'(s)$.
The dual solution is updated to $p' = \prn{0, 0, 1, 0, 0, 0, 0, 0, 0, 1}$ and $q' = \prn{0, 0, 1, 0, 0, 0, 0, 0, 1, 1}$.
Since the corresponding tight coefficient matrix of $A'(s)$ is nonsingular, we stop the algorithm.
The index of the modified DAE remains at two.

Finally, by applying the MS-algorithm to the modified DAE, we obtain an index-1 DAE
\begin{align*}
  \left\{\begin{aligned}
    -i_1 -i_4 +i_5 &= 0, \\
    i_2 + i_3 + i_4 - i_5 &= 0, \\
    v_4 - v_5 &= 0, \\
    z^{[1]}_4 - z^{[1]}_5 &= 0, \\
    -v_1 -v_2 +v_4 &= 0, \\
    v_2 - v_3 &= 0, \\
    R_1 i_1 - v_1 &= 0, \\
    R_2 i_2 - v_2 &= 0, \\
    L \dot{i}_3 - v_3 &= 0, \\
    -i_4 + C z^{[1]}_4 &= 0, \\
    v_5 &= V(t), \\
    z^{[1]}_5 &= \dot{V}(t),
  \end{aligned}\right.
\end{align*}
where $z^{[1]}_4$ and $z^{[1]}_5$ are dummy variables corresponding to $\dot{v}_4$ and $\dot{v}_5$, respectively.

\section{Application to Nonlinear DAEs}
\label{sec:nonlinear_dae}

In this section, we discuss the application of our algorithm to nonlinear DAEs.
The $\sigma \nu$-method~\cite{Chowdhry2004}, which is implemented in Mathematica~\cite{mathematica}, adopts a strategy of treating nonlinear or time-varying terms as independent parameters in the Jacobian matrices of DAEs.
We first describe the $\sigma \nu$-method briefly.

Consider an index-2 nonlinear DAE
\begin{align} \label{def:nonlinear_dae}
  \left\{\begin{alignedat}{4}
    F_1:\:\: & \dot{x}_1 + g(x_2)  &&= f_1(t), \\
    F_2:\:\: & \dot{x}_1 + x_1 + x_3 &&= f_2(t), \\
    F_3:\:\: & \dot{x}_1 \phantom{{}+x_1} + x_3 &&= f_3(t),
  \end{alignedat}\right.
\end{align}
where $\funcdoms{g}{\setR}{\setR}$ is a smooth nonlinear function.
Their method constructs two kinds of Jacobian matrices JD and JV as follows:
\begin{align*}
  \mathrm{JD} = \prn{\pdif{F_i}{\dot{x}_j}}_{i,j} = \begin{pmatrix}
    1 & 0 & 0 \\
    1 & 0 & 0 \\
    1 & 0 & 0
  \end{pmatrix}
  ,
  \quad
  \mathrm{JV} = \prn{\pdif{F_i}{x_j}}_{i,j} = \begin{pmatrix}
    0 & \mathrm{d}g/\mathrm{d}x_2 & 0 \\
    1 & 0 & 1 \\
    0 & 0 & 1
  \end{pmatrix}.
\end{align*}
If JD is nonsingular, the DAE is index zero from the implicit function theorem.
Otherwise, the method performs the Gaussian elimination on JD (and JV simultaneously) to make a row of zeros in JD.
Then the method differentiates the corresponding equation, and checks the nonsingularity of JD again.
The main feature of the $\sigma \nu$-method is to treat nonlinear or time-varying terms as ``independent parameters'' to avoid complicated symbolic manipulations.
The method works according to the rule that arithmetical operations and the differentiation of independent parameters generate new independent parameters.

The $\sigma \nu$-method may fail due to this rule.
For example, let $\alpha_1$ be an independent parameter representing $\mathrm{d}g/\mathrm{d}x_2$ in JV.
By subtracting the first row from the second and third ones, we obtain
\begin{align*}
  \mathrm{JD} = \begin{pmatrix}
    1 & 0 & 0 \\
    0 & 0 & 0 \\
    0 & 0 & 0
  \end{pmatrix}
  ,
  \quad
  \mathrm{JV} = \begin{pmatrix}
    0 & \alpha_1 & 0 \\
    0 & \alpha_2 & 1 \\
    0 & \alpha_3 & 1
  \end{pmatrix},
\end{align*}
where $\alpha_2 = 0 - \alpha_1$ and $\alpha_3 = 0 - \alpha_1$ are newly generated parameters by the rule of arithmetical operations.
We differentiate the second and third rows.
Then JD and JV are
\begin{align*}
  \mathrm{JD} = \begin{pmatrix}
    1 & 0 & 0 \\
    0 & \alpha_2 & 1 \\
    0 & \alpha_3 & 1
  \end{pmatrix}
  ,
  \quad
  \mathrm{JV} = \begin{pmatrix}
    0 & \alpha_1 & 0 \\
    0 & \alpha_4 & 0 \\
    0 & \alpha_5 & 0
  \end{pmatrix},
\end{align*}
where $\alpha_4$ and $\alpha_5$ are parameters corresponding to the derivatives of $\alpha_2$ and $\alpha_3$, respectively.
Although the Jacobi matrix JD is indeed singular due to $\alpha_2 = \alpha_3$, the $\sigma \nu$-method halts at this point as the method regards $\alpha_2$ and $\alpha_3$ as independent.
This failure originates from the elimination of matrices involving the independent parameter $\alpha_1$.
We have confirmed that the implementation in Mathematica actually fails on this DAE.

Our algorithm is applied to the same DAE~\eqref{def:nonlinear_dae} as follows.
Let
\begin{align*}
  A(s) = \begin{pmatrix}
    s   & \alpha &   \\
    s+1 &        & 1 \\
    s   &        & 1 \\
  \end{pmatrix},
\end{align*}
where $\alpha$ is an independent parameter representing $\mathrm{d}g/\mathrm{d}x_2$.
As described in Section~\ref{sec:ms_alg}, the MS-algorithm is applicable to the nonlinear DAE~\eqref{def:nonlinear_dae} if $A(s)$ is upper-tight.
The tight coefficient matrix corresponding to a dual optimal solution $p = \prn{0, 0, 0}$ and $q = \prn{1, 0, 0}$ is
\begin{align*}
  A^\# = \begin{pmatrix}
    1 & \alpha &   \\
    1 &        & 1 \\
    1 &        & 1 \\
  \end{pmatrix},
\end{align*}
which is singular.
Thus we need to modify the matrix.
By the same logic as the discussion in Section~\ref{sec:example_2}, we can regard $A(s)$ as an LM-polynomial matrix $A(s) = \binom{T(s)}{Q(s)}$, where $T(s)$ corresponds to the first row and $Q(s)$ corresponds to the other two ones in $A(s)$.
Then our algorithm modifies $A(s)$ to
\begin{align*}
  A'(s) = \begin{pmatrix}
    s & \alpha &   \\
    1 &        &   \\
    s &        & 1 \\
  \end{pmatrix},
\end{align*}
which is upper-tight (we omit the detail of this modification).
Using an optimal solution $p' = \prn{0, 1, 0}$ and $q' = \prn{1, 0, 0}$ of $\D{A'}$, the MS-algorithm obtains a purely algebraic equation
\begin{align*}
  \left\{\begin{alignedat}{4}
    &z^{[1]}_1 + g(x_2)  &&= f_1(t), \\
    &x_1 &&= f_2(t) - f_3(t), \\
    &z^{[1]}_1 &&= \dot{f}_2(t) - \dot{f}_3(t), \\
    &z^{[1]}_1 + x_3 &&= f_3(t),
  \end{alignedat}\right.
\end{align*}
where $z^{[1]}_1$ is a dummy variable corresponding to $\dot{x}_1$.

This example shows that our algorithm works for a DAE to which the existing index reduction algorithm cannot be applied.
Our algorithm is expected to rarely cause cancellations between nonlinear terms as the algorithm does not perform the row operations involving independent parameters.
Therefore, although the application to nonlinear DAEs remains at the stage of heuristic, we consider that the proposed method can be useful for index reduction of nonlinear DAE.

\section{Conclusion}
\label{sec:conclusion}

In this paper, we have proposed an index reduction algorithm for linear DAEs whose coefficient matrices are mixed matrices.
The proposed method detects numerical cancellations between accurate constants, and transforms a DAE into an equivalent DAE to which the MS-algorithm is applicable.
Our algorithm uses combinatorial algorithms on graphs and matroids, based on the combinatorial relaxation framework.
The algorithm is expected to keep the sparsity of DAEs because the algorithm modifies only rows of accurate constants, and to run much faster in most cases because the algorithm terminates without modifying the DAEs unless the DAEs have unlucky numerical cancellations.
We also have developed a faster algorithm for DAEs whose coefficient matrices are dimensionally consistent.
Our algorithms can be applied to nonlinear DAEs by regarding nonlinear terms as independent parameters.
Numerical experiments on nonlinear DAEs are left for further investigation.

\bibliographystyle{bibstyle}

redpajama_set_name: RedPajamaArXiv


text:
Orthopaedic & Rheumatologic Institute
Orthopaedic & Rheumatologic Institute Menu
Research & Innovations
Home / Institutes & Departments / Orthopaedic & Rheumatologic Institute / Departments & Centers / Pediatric Scoliosis
Our pediatric orthopaedic specialists evaluate and treat cases from the most common to the most complex.
Chat with Appointment Agent
Pediatric Scoliosis patients talk about their experience at Cleveland Clinic.
Scoliosis is an abnormal sideways curvature and rotation of the spine. About 2 to 3 percent of the U.S. population, or 7 million people, have this condition.
Scoliosis may be detected in infancy, childhood or adolescence. While this condition affects both girls and boys, girls are 10 times more likely than boys to need corrective surgery for it.
The thing to pay attention to is whether or not scoliosis progresses. Curves that are larger have a greater chance of progressing than smaller curves. However, even mild curvatures may worsen as a child's spine grows.
Progression can lead to severe consequences in adulthood, including lung and heart problems and potential back pain. That's why it is important to:
Detect scoliosis as early as possible
Monitor its progress
Intervene when necessary
Cleveland Clinic has been providing comprehensive quality care for pediatric patients with acquired and congenital orthopaedic problems for 32 years.
Cleveland Clinic is ranked as one of the nation's top hospitals by U.S. News & World Report.
Cleveland Clinic Children's is also consistently rated among the "Best Children's Hospitals" by U.S. News & World Report.
Our orthopaedic specialists are trained to treat conditions specific to children.
Same day appointment or next business day access is available.
Cleveland Clinic Children's Hospital and other centers have pediatric anesthesiologists on staff. They are adept at managing children's anesthesia, and will provide general anesthesia for spine surgery. Large centers like ours also have a neurological team to monitor spinal cord function minute-by-minute during the entire procedure. After surgery, pediatric anesthesiologists make sure that each child's pain is well-controlled.
View the transcripts from our Online Scoliosis webchats
Download our free pediatric scoliosis treatment guide
There are three types of childhood scoliosis:
Idiopathic scoliosis represents 80 percent of all cases of scoliosis. It usually develops during adolescence, between the ages of 10 and 16. Idiopathic means "of unknown cause" and scientists have not yet unraveled the reasons that scoliosis develops. But because idiopathic scoliosis runs in families, it has a genetic basis.
Congenital scoliosis
Congenital, or infantile, scoliosis is a fairly rare bone abnormality detected at birth. Very often, it accompanies other birth defects, such as heart or kidney problems.
Neuromuscular scoliosis
Neuromuscular scoliosis is caused by abnormalities in the muscles and nerves that support the spine. Thus, patients with muscular dystrophy, spina bifida and cerebral palsy are susceptible. Neuromuscular scoliosis can become quite severe.
We bring together a team of highly trained specialists in order to determine the appropriate intervention and treatment plan for your child. Our pediatric scoliosis team includes:
Pediatric orthopaedic surgeons
Pediatric anesthesiologists
Physical and Occupational therapists
Your child's pediatrician is an integral member of our extended team with whom we consult in developing treatment and, when appropriate, ongoing care plans for your child.
Between our pediatric orthopaedic surgeons' technical expertise and their involvement in ongoing research, we can offer your child the latest, most effective treatments.
Thomas Kuivila, MD
Vice Chair, Orthopaedic and Rheumatologic Institute
Ryan Goodwin, MD
Director for the Center for Pediatric and Adolescent Orthopaedics
David Gurd, MD
Head of Pediatric Spinal Deformity Surgery
To make an appointment with any of our scoliosis specialists, please call the phone number provided below. You may also schedule your appointment online.
Department of Orthopaedic Surgery
Locally 216.444.2606
Toll-free 800.223.2273 ext 42606
Online Second Opinion
We offer second opinions through Cleveland Clinic MyConsult, an effective and secure option for consultation without the time and expense of travel. One of our experts will review your medical records and test results and provide a comprehensive, personalized report, taking care to answer your specific questions.
Learn more about MyConsult Online Medical Second Opinions.
Arthritis & Musculoskeletal Center
The R.J. Fasenmyer Center for Clinical Immunology
Hand, Wrist, Elbow & Shoulder Surgery
Hip Preservation Center
Joint Preservation Center
Musculoskeletal Physical Medicine & Rehabilitation
Orthotic & Prosthetic Services
Osteoporosis and Metabolic Bone Disease Center
Pediatric & Adolescent Orthopaedic Surgery
Rehabilitation & Sports Therapy
Vasculitis Care & Research
View Scoliosis Patient Story

redpajama_set_name: RedPajamaCommonCrawl