Fast Gaussian Elimination for Low Treewidth Matrices

Reducing a matrix to a desired form or decomposing a matrix into a product of matrices with a given structure are basic operations in linear algebra. Problems involving triangular matrices and matrices in row echelon form are often easy to solve, and a large number of tasks may be accomplished by first reducing a general matrix into a matrix of this form. To be precise, a matrix $A=(a_{ij})$ is in row echelon form if all rows consisting entirely of zeros are at the bottom, and the following holds for the remaining rows. If $i$ is a nonzero row, its pivot is the element $a_{ij}\ne 0$ with least $j$. If rows have pivots in columns $j_1$ and $j_2$, respectively, then .

A basic result in linear algebra states that every matrix can be transformed into a matrix in row echelon form by a sequence of three types of row operations, known as elementary row operations. In this paper, we shall refer to two of these types: a type I elementary row operation on a matrix involves interchanging two rows. A type II elementary row operation involves adding a multiple of a row to another row. We say that $A$ and $B$ are row equivalent if $A$ can be transformed into $B$ by a sequence of elementary row operations. The process of reducing a matrix to row echelon form by elementary row operations is known as Gaussian elimination. The standard school algorithm reduces an arbitrary $m\times n$ input matrix $A$ with entries in a field $𝔽$ to row echelon form in $O(mn\min \{m,n\})$ field operations. However, it can be faster if more is known about the structure of $A$. For instance, if $A$ is sparse, its structure can be exploited if we find a pivoting scheme (also known as elimination ordering) that minimizes the fill-in, defined as the set of matrix positions that were initially 0, but became nonzero at some point during the computation. It is worth noting that sparse matrices often occur for some structural reason, such as small treewidth. Gaussian elimination cannot be used efficiently on random sparse matrices with any elimination order, as we would expect such matrices to quickly get dense as the algorithm progresses. We should mention that other strategies can be used to exploit the sparsity of the input matrix, see for instance the work of Peng and Vempala [12].

Much of the early focus has been on the question of characterizing the graphs corresponding to matrices with a perfect elimination order, meaning that there is a pivot strategy for Gaussian elimination without any fill-in. Parter [11] has assigned a graph to a matrix $A$ with an edge $\{i,j\}$ if $a_{ij}$ or $a_{ji}$ are non-zero. He has then shown that there is a perfect elimination order if the graph is a tree. As became customary afterwards, he chooses all pivots in the diagonal and makes the strong assumption that all diagonal elements of the matrix are non-zero and remain so during the computation, usually referred to as not having any accidental cancellation.

Obviously, this assumption simplifies Gaussian elimination enormously. It is justified because it is true for symmetric positive definite matrices, which are the most important special class of matrices for systems of linear equations. However, one would like to see more widely applicable methods.

Rose [16] shows that Parter’s result can be extended to all triangulated graphs (also known as chordal graphs). Rose, Tarjan, and Lueker [18] show that a graph has a perfect elimination order if and only if the graph is triangulated. In this case, they find a perfect elimination order in $O(|V|+|E|)$ time. Rose and Tarjan [17] efficiently find a perfect elimination order in a directed graph if one exists, and show that finding a minimal elimination ordering (minimizing the fill-in) is NP-hard. A directed graph describes the non-zero off-diagonal entries in an $n\times n$ matrix, where Gaussian elimination is done by choosing pivots in the diagonal under the assumption that there will never be a 0 in the diagonal.

The perfect elimination question is more complicated for the case of asymmetric $m\times n$ matrices $A$, where any nonzero pivot is allowed. Any $m\times n$ matrix $A=(a_{ij})$ may be naturally associated with a bipartite graph $G$ with vertex set $\{v_1,\mathrm{\dots },v_m\}\cup \{w_1,\mathrm{\dots },w_n\}$ such that vertices $v_i$ and $w_j$ are adjacent if and only if $a_{ij}\ne 0$. We say that $G$ is the underlying bipartite graph of $A$. This allows us to employ structural decompositions of graph theory to deal with the nonzero entries of $A$ in an efficient way. Different from the symmetric case with pivots restricted to the diagonal, there is no nice characterization known for the bipartite graphs with perfect elimination order when arbitrary pivots are allowed. Nevertheless, a very large class is given by the chordal bipartite graphs defined by Golumbic and Goss [7]. These graphs are not chordal, as obviously only forests are bipartite and chordal. Still, chordal bipartite graphs are characterized similarly to chordal graphs. For chordal graphs, the defining property is that every cycle of length at least 4 has a chord, while chordal bipartite graphs are defined by the requirement that every cycle of length at least 6 has a chord. Golumbic and Goss [7] have shown that chordal bipartite graphs have a perfect elimination order.

They are not the only bipartite graphs with a perfect elimination order, but the only ones where one can always pick any simplicial edge (definition follows) as a pivot [7]. An edge $\{u,v\}$ is simplicial if the union of the neighborhoods of $u$ and $v$ induces a complete bipartite graph. It is easy to see that a bipartite graph $G$ has a perfect elimination order if there exists an ordering of a matching $\{(u_1,v_1),(u_2,v_2),\mathrm{\dots },(u_p,v_p)\}$ in $G$ such that $(u_i,v_i)$ is simplicial in the subgraph $G[V\setminus \{u_1,v_1,u_2,v_2,u_{i-1},v_{i-1}\}]$ for $i\in \{1,\mathrm{\dots },p\}$.

Radhakrishnan, Hunt, and Stearns [13] do Gaussian elimination for symmetric matrices with a given tree decomposition of width $k$ in time $O(k^{2}n)$. As has been common, they implicitly assume that the diagonal is non-zero and there is no accidental cancellation. Naturally, given today’s familiarity with algorithms based on tree decompositions, it would now be an easy exercise to design such an algorithm when pivots can always be selected in the diagonal.

Only in recent years has there finally been some progress towards the difficult problem of efficient Gaussian elimination for matrices of small treewidth without the simplifying assumption that pivots can always be chosen in the diagonal. Fomin et al. [5] have shown a Gaussian elimination algorithm running in time $O(p^{2}n)$ where $p$ is the lesser known tree-partition width parameter or the pathwidth. The tree-partition width of a graph can be arbitrarily large already for graphs of constant treewidth, unless the maximum degree is bounded [19]. However, Fomin et al. [5] still obtained an $O(k^{3}n)$ algorithm for treewidth $k$. We should also refer to the work of Dong, Lee, and Ye [4], whose main result implies an $O(k^{2}n)$ algorithm for the solution of a linear system $Ax=b$ such that the matrix $A\in {ℝ}^{m\times n}$ is full-rank (actually, their result is about the solution of a linear program).

Naturally, the longstanding main goal is an $O(k^{2}n)$ Gaussian elimination algorithm for treewidth $k$ without restrictive assumptions. This would provide a smooth transition to the school algorithm when $k$ approaches $n$. It would also achieve the same running time bound for the general case as for the easy case of non-zero diagonal and no accidental cancellation. It actually seems that the $O(k^{2}n)$ goal has not been explicitly formulated until recently [5], maybe because it had not seemed achievable. However, when this goal has been obtained by simple algorithms under very strong assumptions, the question whether these assumptions (that convenient entries are nonzero and thus can be used as pivots) are necessary must have been obvious.

We make substantial progress in this direction by devising the algorithm Fast Gaussian Elimination. Given an $m\times n$ matrix $A$ over a field $𝔽$ with underlying bipartite graph $G$, and given a tree decomposition $𝒯$ of $G$ of width $k$ and size $O(k(m+n))$, it produces a matrix $U$ in row echelon form that is row equivalent to $A$. The entries of the columns of $U$ that contain a pivot (and therefore generate a space of dimension $rank(A)$) are fully computed in time $O(k^{2}n)$. The entries of the remaining columns of $U$ are described implicitly as linear combinations of at most $k$ other columns, and may be computed explicitly in time $O(km(n-rank(A))))$. In particular, this algorithm allows us to compute the rank and the determinant of $A$ in time $O(k^2(m+n))$. It also allows us to decide in time $O(k^2(m+n))$ whether the linear system $Ax=b$ has a solution, and to compute a solution of the linear system in case it exists. Our main result is stated more formally below.

We actually also compute a permutation matrix $P$ of order $m$ such that $U$ is the row echelon form that would be obtained from $PAQ$ by the most standard algorithm (adding multiples of earlier rows to later rows).

Our algorithm is not based on the $O(k^{3}n)$ algorithm of Fomin et al. [5], as it uses standard tree decompositions in a natural way, and does not rely on the somewhat obscure notion of tree-partition width. Even though our algorithm removes the decades old severe restriction to pivoting in the main diagonal, the main idea of our Gaussian elimination algorithm is quite simple. The tree decomposition provides an elimination order for rows (equations) and columns (variables). This order sequentially hands out permissions to rows and columns to be eliminated. However, when a row $i$ receives a green light to be processed, it is only processed immediately if it contains a variable that has already received a green light. Otherwise, this row is put into a buffer to be delayed until one of its variables $x_j$ also has a green light. Then $a_{ij}$ is chosen as a pivot. Likewise, when a variable receives a green light then, if needed, its column is put into a buffer until some row containing this variable has a green light.

As was the case in Fomin et al. [5] we use a nice tree decomposition of the underlying bipartite graph $G$ as defined by Kloks [9]. However, to deal with $0$’s at intended pivot locations, we create temporary buffers for rows and columns that allow us to delay the selection of pivots as discussed above. As is typical for algorithms based on tree decompositions, row operations are actually performed on small submatrices of the input matrices. For simplicity, the algorithm does not manipulate the large given matrix. Instead, these small submatrices, called boxes, are maintained separately. The small submatrices in this paper include additional information to keep track of the rows and columns that lie in the temporary buffer, which allows us to make headway towards the output matrix in row echelon form even when accidental cancellations occur. It is crucial to perform so-called bookkeeping operations that ensure that the temporary buffers remain small during the entire application of the buffer.

The paper is organized as follows. The algorithm is described in Section 2, where we also discuss its correctness. The running time is analyzed in Section 3. We should mention that, to achieve this running time, we need to keep track of the field operations in a global way, rather than adding the contributions of the worst-case scenarios at each step. Together, these two sections establish Theorem 1. We provide several consequences of our result in Section 4.

We now describe the algorithm Fast Gaussian Elimination. We start with the definition of a tree decomposition, which we now state. Let $G=(V,E)$ be a $\nu$-vertex graph, with the standard assumption that $V=[\nu ]=\{1,\mathrm{\dots },\nu \}$. A tree decomposition of a graph $G$ is a tree $𝒯$ with nodes $\{1,\mathrm{\dots },r\}$, where each node $t$ is associated with a bag $B_t\subseteq V$, satisfying the following properties: (1) ${\bigcup }_{t=1}^r{B}_t=V$; (2) For every edge $\{v,w\}\in E$, there exists $B_t$ containing $v$ and $w$; (3) For any $v\in V$, the subgraph of $𝒯$ induced by the nodes whose bags contain $v$ is connected. The width of the tree decomposition $𝒯$ is defined as ${\max }_t\left(|B_t|-1\right)$ and the treewidth $\mathrm{tw}(G)$ of graph $G$ is the smallest $k$ such that $G$ has a tree decomposition of width $k$. Tree decompositions were popularized by the seminal work of Robertson and Seymour [14, 15], but other definitions that are similar or even equivalent have appeared in earlier work, see [1, 8]. As mentioned in the introduction, we use the concept of nice tree decomposition introduced by Kloks [9], which is a rooted tree decomposition $𝒯$ of a graph $G$ such that all nodes are of one of the following types:

1. (a)

   (Leaf) The node $t$ is a leaf of $𝒯$;

2. (b)

   (Introduce) The node $t$ introduces vertex $v$, that is, it has a single child $s$, $v\notin B_s$ and $B_t=B_s\cup \{v\}$.

3. (c)

   (Forget) The node $t$ forgets vertex $v$, that is, $t$ has a single child $s$, $v\notin B_t$ and $B_s=B_t\cup \{v\}$;

4. (d)

   (Join) The node $t$ is a join, that is, it has two children $s$ and $s^{\prime }$, and $B_t=B_s=B_{s^{\prime }}$.

We further assume that the bag associated with the root is empty. This is easy to accomplish, as any tree decomposition $𝒯$ that satisfies (a)-(d), but does not satisfy this additional requirement, may be turned into a nice tree decomposition with empty root bag by appending a path of length at most $k+1$ to the root of $𝒯$, where all nodes have type forget. Having a root with an empty bag ensures that each vertex of $G=(V,E)$ is associated with exactly one forget node that forgets it. Moreover, it ensures that, for any $v\in V$, the node of $𝒯$ that is closest to the root among all $t$ such that $v\in B_t$ is the child of the node that forgets $v$. Despite the additional structure, a nice tree decomposition may be efficiently derived from an arbitrary tree decomposition. More precisely, Kloks [9, Lemma 13.1.2] has shown that, if $G$ is a graph of order $\nu$ and we are given an arbitrary tree decomposition of $G$ with width $k$ and $r$ nodes, it is possible to turn it into a nice tree decomposition of $G$ with at most $4\nu$ nodes and width at most $k$ in time $O(k(\nu +r))$. Regarding the computation of a tree decomposition, Bodlaender’s theorem [2] provides an optimal tree decomposition in linear time when the treewidth $k$ is a constant. Bodlaender et al. [3] and Korhonen [10] compute faster constant factor approximations in linear time. Thus, for bounded $k$, there is no need to assume a tree decomposition is given, since both the tree decomposition and the Gaussian elimination are then computed in time $O(n)$. However, the constants hidden in the $O$-notation grow quickly with $k$.

We are now ready to describe the algorithm. Let $A$ be an $m\times n$ matrix with entries in a field $𝔽$ and let $G=(V,E)$ be the underlying bipartite graph with vertex set $V=\{v_1,\mathrm{\dots },v_m\}\cup \{w_1,\mathrm{\dots },w_n\}$ associated with it. We wish to find an $m\times m$ permutation matrix $P$, an $n\times n$ permutation matrix $Q$ and an $m\times n$ matrix $U$ in row echelon form such that $PAQ$ is row equivalent to $U$ in a direct way. We shall operate on a nice tree decomposition $𝒯$ of $G$ with node set $[r]$ and width $k$ rather than on the matrix $A$ itself. The algorithm Fast Gaussian Elimination works bottom-up on the rooted tree $𝒯$, that is, it only processes a node $t$ after its children have been processed. Each node $t$ except the root produces a data structure known as a box and transmits it to its parent. A box is a matrix $A_t$ with $O(k)$ rows and columns. It may be recorded as a triple of matrices $(M_t,N_t,J_t)$ whose rows and columns are labeled by distinct elements of $\{v_1,\mathrm{\dots },v_m\}$ and $\{w_1,\mathrm{\dots },w_n\}$, respectively, associating rows and columns of the boxes with rows and columns of the original matrix. We may view it as

where $J_t$ is an $r_t\times c_t$ matrix whose rows and columns are labeled by the elements in $B_t$, the bag associated with node $t$. Here $r_t$ is the number of row vertices (i.e., vertices $v_i$) in $B_t$ and $c_t$ is the number of column vertices (i.e., vertices $w_j$) in $B_t$. The matrices $M_t$ and $N_t$ have dimensions $m_t\times c_t$ and $r_t\times n_t$, respectively, where and . Observe that some of these matrices can be degenerated, in the sense that $m_t=0$, $n_t=0$, $c_t=0$ or $r_t=0$. In a case where $m_t=0$, $n_t=2$, $c_t=2$ and $r_t=3$, for instance, the box $A_t$ is $3\times 4$ matrix, $M_t$ is a matrix with no rows and two empty columns, and both $N_t$ and $J_t$ are $3\times 2$ matrices. We should note that a similar data structure, also called a box, has been exploited in the context of computing a diagonal matrix that is congruent to a given symmetric matrix using a tree decomposition of its underlying graph [6].

It will also be useful to assign a canonical ordering to the rows and columns of the input matrix. We assume that they are ranked according to the order in which their Forget nodes appear in some topological ordering (e.g., obtained by a post-order traversal of $𝒯$), that is, $v_1$ is the first vertex that is forgotten among all vertices that represent rows, $v_2$ is the second, and so on. The same applies to $w_1,\mathrm{\dots },w_n$.

We start with an intuition about the meaning of the boxes and about how the algorithm works. Recall that the nodes of the nice tree decomposition $𝒯$ are ordered bottom-up as $1\mathrm{\dots },r$, where $r$ is the root and $B_r=\mathrm{\varnothing }$. At each node $t\in V(𝒯)$, the algorithm initializes box $A_t$ (if $t$ is a leaf) or produces a new box $A_t$ based on its bag $B_t$ and on the boxes transmitted by its children. While producing the box $A_t$, the algorithm may access entries $a_{ij}$ of the input matrix such that both $i$ and $j$ are associated with vertices in $B_t$. The box $A_t$ is then transmitted to its parent, except if $t$ is the root of $𝒯$. We call this processing node $t$.

It is also useful to keep in mind that the algorithm performs two types of elementary row operations, which we call authentic operations and bookkeeping operations. The former are operations that can be performed in the full matrix $A$ in order to reach row echelon form. The latter are operations that are only performed on the boxes as a way to limit their sizes, but which do not have full matrix counterparts. The effect of processing nodes may be viewed as a sequence of $m\times n$ matrices

where $U$ is in row echelon form and ${\tilde{A}}_t$ is obtained from ${\tilde{A}}_{t-1}$ by performing the authentic elementary row operations that have been performed while node $t$ is processed. We should mention that this sequence of row equivalent matrices is a simplification of what actually happens, but that it could be achieved if an oracle gave us the proper ordering of rows and columns of $A$. In reality, the proper orderings are also computed by the algorithm while the nodes are processed, leading to the permutation matrices $P$ and $Q$ mentioned in the statement of Theorem 1.

To achieve the complexity stated in Theorem 1, row operations cannot be actually performed on $m\times n$ matrices, but are instead performed only in the (smaller) boxes. It is useful to think that any row or column of the input matrix $A$, say row $i$, is initially untouched and changes to a regular row when the vertex $v_i$ associated with it appears in the bag of a node of the branch¹¹1We use “branch” to refer to the subtree induced by all descendants of a node. of the tree decomposition rooted at $t$. When $v_i$ is removed from the bag (at the node that forgets it in the nice tree decomposition), either row $i$ is assigned a pivot in some column $j$, in which case row $i$ is classified as processed (column $j$ is also classified as processed), or row $i$ becomes a buffer row, meaning that it has been given the green light for a pivot, but that the decision of assigning a pivot has been deferred to a later step because no column has received a green light. In the process of keeping the size of the box bounded, the algorithm may decide that it is impossible or unnecessary to assign a pivot to a buffer row or to a buffer column, in which case the status of the row or column also changes to processed. We observe that the classification of each row and column is defined independently for disjoint branches of the tree decomposition, and an important feature of the algorithm is that these independent classifications can be made compatible when two branches are merged (which happens at join nodes).

This finally gets us to the meaning of a box $A_t=(M_t,N_t,J_t)$ at node $t$. The rows and columns of the matrix $J_t$ are precisely the regular rows and columns of the branch ${𝒯}_t$ of $𝒯$ rooted at node $t$. In other words, these are the rows and columns corresponding to the vertices in the bag $B_t$. Moreover, each entry $ij$ of $J_t$ records the net change to the original entry $a_{ij}$ due to row operations performed while processing nodes in ${𝒯}_t$. The matrix $M_t$ records the entries $ij$ of the matrix ${\tilde{A}}_t$ defined in (2) in the case where $i$ is a buffer row (with respect to ${𝒯}_t$) and $j$ is a regular column. Analogously, $N_t$ records the entries $ij$ of ${\tilde{A}}_t$ in the case where $i$ is a regular row and $j$ is a buffer column. This reveals another important feature of the algorithm: The entries of the full $m\times n$ matrix associated with regular rows and regular columns can easily be modified in parallel by computations in disjoint branches, as the algorithm does not capture their exact values while processing a single branch. The boxes just store the net changes to the values due to operations in this branch. However, this does not happen for buffer rows and buffer columns, for which exact values are recorded in the box. The top left corner of $A_t$ in (1) shows that, if $i$ is a buffer row and $j$ is a buffer column, the entry $ij$ in ${\tilde{A}}_t$ is equal to zero. This is consistent with the idea that the algorithm selects a new pivot whenever there are a buffer row $i$ and a buffer column $j$ for which the element $ij$ is nonzero, as both the row and the column had already been given the green light as possible rows and columns for a pivot. So, if row $i$ and column $j$ are in the buffer at the end of step $t$, then the entry $ij$ must be 0.

While producing the boxes, the algorithm may select pivots for the output matrix $U$, it may identify zero rows that cannot contain a pivot and it may identify columns that are linear combinations²²2Actually, they are not necessarily linear combinations of the other columns with respect to the entire matrix, but they are linear combinations with respect to a submatrix as depicted in Figure 1. of other columns, and therefore cannot contain pivots together with the other columns. The status of these rows and columns changes to processed, which means that information about them is recorded by global variables, while they are removed from the box.

After this description, we come back to a matrix ${\tilde{A}}_t$ in (2), depicted in Figure 1³³3The authors thank Elizandro Max Borba for the figure.. The top left entries lie in processed rows and processed columns that contain pivots. These entries are not modified later in the algorithm. The bottom rows are (processed) zero rows that have been produced while limiting the size of the set of buffer rows (it is indeed impossible to find a pivot in such a row), while the rightmost columns are (processed) columns that have been identified as linear combinations of other columns (except for possible entries in rows whose pivot had already been computed when this column was processed). As a consequence, it is unnecessary to assign a pivot to them. Moreover, while the actual values in the rightmost columns may be modified in later steps of the algorithm, we need not keep track of them, as each such column may be computed as a linear combination of (a known set of) $O(k)$ columns with smaller index at the end of the algorithm. The central submatrix represents entries in the buffer, and in regular or untouched rows and columns. Some of them lie in the box where computations are actually performed in a given step.

Before giving a formal description of every step of the algorithm, we briefly describe what is done at a node of each type. A node of type leaf merely initializes a box with regular rows or columns, while a node of type introduce simply adds a new regular row or column to the box transmitted by its child. Nodes of type join combine the information that is passed on by two different branches of the decomposition, which have been processed independently of each other. Let $A_j$ and $A_{\mathrm{\ell }}$ denote the boxes produces by the children of $t$. A new matrix $M_t$ is produced by stacking the matrix $M_j$ on top of $M_{\mathrm{\ell }}$, and a new matrix $N_t$ is produced by juxtaposing $N_j$ to the left of $N_{\mathrm{\ell }}$. Such matrices $M_t$ and $N_t$ may become too large, so that row operations are performed to remove some rows and/or columns from the buffer. Producing $J_t$ also requires merging $J_j$ and $J_{\mathrm{\ell }}$ , as the regular rows and columns are the same. Being able to easily merge the matrices transmitted by the children is the reason why the rows and columns are considered with the canonical ordering mentioned above.

Nodes of type forget concentrate most of the action. Recall that each row $i$ of $A$ is associated with a vertex $v_i$ of $G$, which in turn is associated with the unique Forget node of $𝒯$ that forgets $v_i$. Analogously, each column $j$ of $A$ is associated with a vertex $w_j$ of $G$ and the unique Forget node of $𝒯$ that forgets $w_j$. Assume that the algorithm is processing a Forget node $t$, which means that it received a box $A_s$ of its child $s$. Further assume that node $t$ forgets a vertex $v_i$ corresponding to row $i$. This means that the algorithm singles out this row as a candidate for having a pivot. In the box $A_s$, it checks whether there is a column $j$ that has already been singled out for a possible pivot, but which has no pivot so far. These are precisely the columns of $N_s$. If such a column exists, it checks whether the entry $ij$ in $N_s$ nonzero. If this is also satisfied, the algorithm performs operations to create a pivot in row $i$ and column $j$, and they are both removed from the box. If one of the conditions fails, row $i$ is removed from $J_s$ and added to $M_t$, which contains rows singled out for a pivot, but which have no pivot so far. Furthermore, if the condition on the maximum number of rows of $M_t$ fails after the addition of this new row, the algorithm acts to remove some rows from the buffer.

Analogously, if node $t$ forgets vertex $w_j$ corresponding to column $j$, it first checks whether there are rows that have been singled out for pivots, but have not been used, namely the rows of $M_s$. It then looks for one such row $i$ for which the entry $ij$ is nonzero. If such a row exists, the algorithm performs operations to create a pivot in row $i$ and column $j$, and they are both removed from the box. If one of the conditions fails, column $j$ is removed from $J_s$ and added to $N_t$, which contains columns singled out for a pivot, but which have no pivot so far. If the condition on the number of columns of $N_t$ fails after the addition of this new column, the algorithm acts to remove some columns from the buffer.

A high level description of the algorithm may be found in Figure 2. Note that, in this descriptiom, we have five global variables: the number $\rho$ of pivots already computed, two vectors $𝐮$ and $𝐩$ of length $m$ and two vectors $𝐰$ and $𝐪$ of length $n$. Initially, these are zero vectors. At the end of the algorithm, vector $𝐮$ allows us reconstruct the rows of the output matrix $U$, while vectors $𝐩$ and $𝐪$ encode the permutations that generate $P$ and $Q$, respectively. As we discussed, the algorithm may also find columns of $A$ that are linear combinations of other columns and remove them from boxes. When this happens, the linear combination that produces it is recorded in the vector $𝐰$, so that the values of entries of $U$ in these columns can also be recovered at the end of the algorithm.

In the remainder of this section, we shall describe each operation in detail. To simplify our discussion, for a node $t$, we say that buffer rows and buffer columns have type 1, while regular rows and regular columns have type 2.

In this section, we find an upper bound on the number of operation required to produce matrices $P$, $Q$ and $U$ as in the statement of Theorem 1 using algorithm Fast Gaussian Elimination. If the algorithm is given a tree decomposition $𝒯$ with $r$ nodes and width $k$ of the underlying graph $G$ of matrix $M$, it first computes a nice tree decomposition ${𝒯}^{\prime }$ of the same width $k$ with fewer than $4(m+n)$ nodes in time $O(k(r+m+n))$ as discussed above. As it turns out, the number of nodes of each type in ${𝒯}^{\prime }$ is at most $m+n$.

We consider the number of operations performed at each type of node. Leaf Box initializes a matrix in $r_t\times c_t=O(k^2)$ trivial steps. Introduce Box uses $O(k)$ steps to create a row or column filled with zeros.

For Join Box and Forget Box, the main cost comes from row operations. As each row has size at most $n_t+c_t\le 2r_t+c_t\le 2(r_t+c_t)\le 2k+2$, each such row operation requires $O(k)$ additions and multiplications. Regarding Forget Box, when vertex $v_i$ associated with a row $i$ or vertex $w_j$ associated with a column $j$ is forgotten, either a new type 1 row or column is created, or a new pivot is found. If the latter happens, the algorithm performs at most $m_t+r_t-1=O(k)$ row operations, where each row has size at most $O(k)$, leading to $O(k^2)$ field operations. Additionally, the algorithm may perform row operations to keep $M_t$ or $N_t$ of the right size, just as in Join Box (see below). If a new type 1 row or column is created, the algorithm just inserts the new row in $M_t$ or the new column in $N_t$, possibly performing row operations to keep $M_t$ or $N_t$ of the right size, just as in Join Box.

It remains to account for the row operations performed to keep matrices $M_t$ and $N_t$ of the right size. Recall that, every time row operations of this type are necessary to produce $M_t$, we start with a matrix of dimension $m^{\ast }\times c_t$, where $m^{\ast }\ge 2c_t$ and we create at least $m^{\ast }-c_t\ge c_t$ zero rows. These rows become processed. For $N_t$, we start with a matrix of dimension $r_t\times n^{\ast }$, where $n^{\ast }\ge 2r_t$ and we identify at least $n^{\ast }-r_t\ge r_t$ columns that are linear combinations of other columns. Again, the columns become processed. We shall keep track of the overall number of field operations performed at these steps by “charging” the rows and columns processed at this step, that is, if $f$ field operations are performed and $\mathrm{\ell }$ rows are processed, each of the $\mathrm{\ell }$ rows is charged with $f/\mathrm{\ell }$ operations. Because of the dimensions of the matrices, getting them to row echelon form takes $O(m_{}^{\ast }{}_{}{}^2{c}_t)$ and $O(n^{\ast }{r}_t^2)$ field operations, respectively. In particular, each row can be charged at most $O(m^{\ast }{c}_t)=O(k^2)$ and each column can be charged at most $O(r_t^2)=O(k^2)$. Since any row or column may be processed at most once in this way, the overall number of field operations that the algorithm performs to keep matrices $M_t$ or $N_t$ of the right size is $O((m+n)k^2))$.

Overall, the number of operations taken by the algorithm to produce the vectors $𝐩$ and $𝐪$, in addition to vectors $𝐮(i)$ and $𝐰(j)$ for every $i\le \rho =rank(A)$ and every $j>\rho$ is $O(k(|𝒯|+m+n)+k^2(m+n))$.

We now consider the decomposition $U_t=P_t{\tilde{A}}_t{Q}_t$. The matrices $P$ and $Q$ are defined using $𝐩$ and $𝐪$ with no additional field operations. Indeed, $P=(p_{ij})$ is the permutation matrix such that $p_{𝐩(i)i}=1$ for $1\le i\le m$, while $Q=(q_{ij})$ is the permutation matrix such that $q_{𝐪(i)i}=1$ for $1\le i\le n$.

For $U$, the entries may be computed in two ways. They are either obtained from a vector $𝐮(i)$ with no additional computation, or they are obtained as a linear combination of $O(k)$ earlier values using coefficients recorded in $𝐰(j)$, which requires $O(k)$ sums and products. The number of operations performed to obtain entries of the second type is at most

This concludes the analysis.

We conclude the paper with the statement of results that are immediate consequences of Theorem 1.