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{equation} \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} \end{equation}

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.

Thank you in advance.

  • 1
    $\begingroup$ 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. $\endgroup$ – LSpice Apr 13 at 22:01
  • $\begingroup$ 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. $\endgroup$ – avril_14th Apr 13 at 22:37
  • $\begingroup$ 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? $\endgroup$ – Federico Poloni Apr 15 at 13:19
  • $\begingroup$ 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. $\endgroup$ – avril_14th Apr 15 at 20:22
  • $\begingroup$ 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. $\endgroup$ – 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 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.

  • $\begingroup$ 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). $\endgroup$ – avril_14th Apr 14 at 17:09
  • $\begingroup$ 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). $\endgroup$ – avril_14th Apr 14 at 17:16
  • $\begingroup$ 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. $\endgroup$ – Federico Poloni Apr 14 at 17:38
  • $\begingroup$ 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. $\endgroup$ – Federico Poloni Apr 14 at 17:42
  • $\begingroup$ 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. $\endgroup$ – avril_14th Apr 14 at 17:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.