# The “best way” to order unknowns in linear systems

Start with a linear system of the form $$\begin{equation*} Ax + Bt + C = 0, \end{equation*}$$ where $$x = (x_1, \dots, x_n) \in \mathbb R^n$$ is the vector of unknowns, $$t \in \mathbb R^m$$ is a vector of parameters, $$A \in GL(n, \mathbb R)$$, $$B \in \mathcal M_{n,m}(\mathbb R)$$ and $$C \in \mathbb R^n$$. Suppose you can apply Gauss reduction without pivoting. Then you end up with a system of the form $$\begin{equation*} Ux + Pt + Q = 0, \end{equation*}$$ where $$U$$ is a $$n \times n$$ unitriangular matrix, $$P \in \mathcal M_{n,m}(\mathbb R)$$ and $$Q \in \mathbb R^n$$. This allows to compute $$x$$ recursively by: $$$$\begin{cases} x_n = -(p_{n1} t_1 + \dots + p_{nm} t_m + q_n) \\ x_{n-1} = -(u_{n-1,n} x_n + p_{n-1,1} t_1 + \dots + p_{n-1,m} t_m + q_{n-1}) \\ \dots \\ x_1 = -(u_{12} x_2 + \dots + u_{1n} x_n + p_{11} t_1 + \dots + p_{1m} t_m + q_1) \\ \end{cases}$$$$

Of course, you can do this for every permutation of $$(x_1, \dots, x_n)$$, and also for every permutation of rows that does not give you zeroes on the diagonal at some point of Gauss reduction. Now my questions are:

1. Is it possible to know in advance (i.e. without doing the whole computation) which permutation of rows does not make Gauss reduction fail because of zeroes on the diagonal? I think this is very classical but I did not find any reference.

2. Is it possible to know in advance which permutation of rows and columns gives you the triple $$(U, P, Q)$$ with the greatest number of zeroes? This would give me the best algorithm that computes $$x$$ from a given $$t$$.

3. Does it help if you make some assumptions on the form of $$A$$? In my situation $$A$$ is sparse, often it is also symmetric and weakly diagonally dominant.

• What purpose does $t$ serve? That is, are we solving for $t$ (in which case it seems that it's really part of the vector of variables) or regarding $t$ as fixed (in which case it seems that $B t + C$ is just a complicated way of naming the inhomogeneity)? It seems like you're asking for a quick way to determine to which Bruhat coset $B w B$ the matrix $A$ belongs, but I'm not sure. – LSpice Apr 13 at 22:01
• I edited the question because it was partially incorrect and I also gave more details. The vector $t$ is a parameter and I'm searching for the best algorithm that computes $x$ for a given $t$. If you want, you can think of $t$ as the "input" and $x$ as the "output" of the system. – avril_14th Apr 13 at 22:37
• Do you have many different choices of $t$ that you need to compute the solution for? How many? Could you give an idea of the size of the numbers involved? – Federico Poloni Apr 15 at 13:19
• Essentially, my aim is writing very fast code to solve discrete dynamical systems: the situation I described is what happens in the linear case. The vector $t$ depends on the previous step and in principle it could be any vector, but I would say that in my concrete problem $|t_i| \leq 10^5$ for every $i$ and $10^{-5} \leq |a_{ij}| \leq 10^5$ whenever $a_{ij} \neq 0$. Using some heuristics I can already write stable code for all situations I'm interested in, now I'm trying to improve efficiency. – avril_14th Apr 15 at 20:22
• I meant matrix sizes and sparsities (and how many different values of $t_i$ there are), not the magnitude of their entries. Sorry for the sloppy wording from my side. – Federico Poloni Apr 16 at 6:44

(1): the classical condition is that all leading principal submatrices of $$A$$ (apart from $$A$$ itself) are nonsingular; see e.g. Theorem 5.2.3.4 here. This applies because what you are computing is an LU factorization of $$A$$ (I see you don't use this term; if it is not familiar to you, I suggest you start from reading about it, because it is a very useful concept for this kind of problems.) This condition might settle some special cases, but in general as far as I know there is no better way to tell if the leading principal matrices are singular than computing the factorization and seeing if it succeeds.

(2) Optimally, no. Even with $$P=Q=0$$ and $$A$$ positive-definite, computing the minimum fill-in is NP-complete. There are various heuristics though, e.g. this, because it is a common problem in numerical linear algebra. Basically every sparse LU package contains such a heuristic, so normally as an end-user you don't have to worry about that.

(3) Not really, as stated in (2). Usually you worry about fill-in only for sparse matrices to begin with, and even for positive definite matrices (which are diagonally dominant up to a scaling that does not affect sparsity) the NP-completeness result holds.

Addendum: I don't think that this is going to give you the best algorithm to compute $$x$$ given a $$t$$, as you claim in (2). Rather, I would, given a $$t$$, (a) evaluate $$Bt+C$$ to a parameterless matrix (b) apply to it the inverses of the elimination matrix $$L$$ and $$U$$ obtained by Gaussian elimination. In this way, you can disregard sparsity in $$P,Q$$ and focus on reducing fill-in in $$U$$ only.

• 1) I know about the classical condition, but I was searching for some easier criterion that allows you to know in advance. I agree with you that this may be not possible. 2) Thank you, this is the answer I was searching for. Can you tell me more about available heuristics? I don't need to solve specific instances of the problem, but rather to write code that works for all possible values of $t$, and the entries of my matrices are parametric (so I cannot simply run a numerical routine on $A$ once for all). – avril_14th Apr 14 at 17:09
• Regarding the addendum, if I want to reduce the total number of operations to a minimum I cannot ignore the structure of $B$ and $C$, as computing $Bt + C$ may be expensive. After thinking more about the problem, I found a much more elegant formulation: given an affine map $x: \mathbb R^m \to \mathbb R^n$, $t \mapsto At + B$, which is the best algorithm to evaluate it? Of course a non-recursive algorithm has in general quadratic complexity, but you can improve it by allowing recursions on $x$. I think this is a very hard problem (it includes reduction from DFT to FFT for example). – avril_14th Apr 14 at 17:16
• I am not an expert on sparse reordering methods; unfortunately it gets quite technical. Some first pointers are arxiv.org/abs/1907.05309 and the books "Direct Methods for Sparse Linear Systems" by Tim Davis and "Direct Methods for Sparse Matrices" by Duff et al. – Federico Poloni Apr 14 at 17:38
• Does $A$ depend on $t$ as well? It didn't seem so from your statement. If it does, unfortunately I don't think it's going to be simple to write code that works robustly for all $t$; it may be better to recompute the factorization from scratch for each value of $t$. Similarly, if $B$ and $C$ are not sparse and do not have a very strong specific structure such as your example of FFT matrices, I don't think there is much hope to spend less than $O(n^2)$ to compute $Bt+C$. Recursions will work only for very specific problems. – Federico Poloni Apr 14 at 17:42
• No, $A$ does not depend on $t$. As you said, if my matrices don't have a specific structure there is no hope to reduce asymptotic complexity, but in my "concrete" situation even saving a few operations might be useful. On the other hand, usually $n \leq 12$, so in the majority of cases brute force will work. I was wondering if there is at least some criterion to reduce possible solutions or some greedy algorithm. – avril_14th Apr 14 at 17:57