13

In general, there is no shortcut other than completely re-factoring the matrix. There have been a few similar questions on this SE that cover the topic in more depth than I can: Can diagonal plus fixed symmetric linear systems be solved in quadratic time after precomputation? LU Decom of PSD Matrix + Diagonal Matrix Perturbation of Cholesky decomposition ...


7

I think you are on the right way. If you correct your errors, it will look very similar to http://www.math.toronto.edu/mpugh/Teaching/Mat1062/notes2.pdf.


7

Yes. The block Lanczos algorithm http://www.netlib.org/utk/people/JackDongarra/etemplates/node250.html produces a block triangular matrix where you control the block size, hence the bandwidth. Certainly, one can argue that a block tridiagonal matrix is not a "proper" banded matrix as there regular patches of certain zeros within the band. If you want ...


6

If the only non-zero entries of $A_{ij}$ have $j$ in $\{i - 1, i, i + 1\}$, then $A$ is a banded matrix with bandwidth 1. More generally, you can talk about matrices of bandwidth $k$ where $k$ is any integer. For example, if you were using a higher-order finite difference discretization that used more points to calculate a derivative, you'd get higher-...


5

If you are trying to calculate all 50k eigenvalues of these matrices, that is going to take you a "quite good amount of time" no matter which algorithm you choose. I'm doubtful you can do substantially better than LAPACK. Lanczos (or generally Arnoldi) algorithms are particularly good at calculating a hand full of eigenvalues at either end of the spectrum (...


5

Nothing is wrong here. The Cuthill-McKee algorithm is a greedy algorithm, and doesn't depend too much on the order of A on input. The reverse Cuthill-McKee algorithm is often used to produce nice orders for skyline solvers, and the skyline of the reordered matrix looks indeed quite reasonable. (The bandwidths of Cuthill-McKee and reverse Cuthill-McKee is the ...


4

Well, other than the usual "don't invert your matrices unless you need the inverse itself" you can still use the banded routines ?gbtrf and then use ?gbtrs with the right hand side being an identity matrix.


3

Computing an approximation of $P_1=e^Q$ is typically done by some kind of polynomial expansion of the matrix exponential. One doesn't take the Taylor expansion (it converges too slowly) but for the sake of the argument, let me assume that you do and that you approximate $$ P_1 = e^Q \approx \sum_{k=0}^N \frac1{k!} Q^k = I+Q+\frac 1{2!} Q^2 + \frac 1{3!}Q^...


2

There are many matrix reordering algorithms, and they seek a permutation array that maps from the original coordinates to the reordered coordinates. Behrisch et al. (2016) gave a nice review of different reordering algorithms, with an analysis of their usage, complexity, objective and applications. The objective is not always to minimize matrix bandwidth, ...


2

A software package mentioned in Moler and van Loan's "twenty-five years after" portion of their paper (Sec. 13, see link in Christian Clason's comment) is Roger Sidje's Expokit, with an accompanying paper that explains the approach. For large sparse matrices a Taylor-like polynomial is applied to Krylov subspaces, obtaining the advantages of matrix-vector ...


2

The special structure of your matrix is easily exploited by a custom made $QR$ factorization based on Givens rotations. This leads to a method of $O(n)$ complexity. Below is a Octave/Matlab code, which solves $Mx=rhs$ for cyclic bidiagonal matrices $M$ of size $n \ge 2$ . % Create a cyclic bidiagonal testmatrix a=[2 -3 -1 5 3]'; b=[7 2 -4 -1 6]'; rhs=[-1 -3 ...


2

If you can solve the linear system with a direct solver, then that's exactly what you should be doing. Multigrid is a method that can be used if you don't have the time or memory resources to use a direct solver (because direct solvers have a complexity that grows faster than $O(N)$ with the size of a linear system). If you use a direct solver as a sub-step ...


2

@gohokies has already given the correct answer in a comment, but just for more context: Matlab backslash calls the UMFPACK (now SuiteSparse) solver for sparse linear systems. The default ordering used by UMFPACK is indeed the Approximate Minimum Degree (AMD) method, or a variation thereof.


1

1) Is there a mathematical trick to simplify the above matrix equations? As in, only having to do one inverse operation instead of two. Yes: Schur complement formulation. Your system is equivalent to the larger one $$ \begin{bmatrix} 0 & B\\ B^T & A^T \end{bmatrix} \begin{bmatrix} y\\ z \end{bmatrix} = \begin{bmatrix} -b\\0 \end{bmatrix} $$ with $b=...


1

I would assume, that your $N\times N$ system has the following form: $$ \underbrace{\begin{pmatrix} b_1 & c_1 & & & & 0 \\ a_1 & b_2 & c_2 & & & \\ s_1 & a_2 & b_3 & c_3 & & \\ & s_2 & a_3 & b_4 & \ddots & \\ &...


1

The equation: $$\partial_t f + \partial_xF = 0$$ with $f_{x=0} = f_0(t)$ (which is known) and appropiate initial conditions can be discretised as follows: $$\partial_t\vec{f}+D\vec{F}=\vec{b}\tag{*}$$ Where $\vec{f} = (f_1,f_2,...,f_{N+1})^T$ and $D$ is the discrete operator which replaces $\partial_x$, and is equal to: $$D = \frac{1}{h}\left[\begin{matrix}...


1

From your comment "The plot emphasises expectedly the connections of node 1", I guess that maybe the idea is showing that node 1 is connected not only to a set of nodes that are "one close to the other", but its interconnections are spread somewhat evenly across the graph. Reverse Cuthill-McKee in theory reorders the nodes so that clusters are mapped into ...


1

The most efficient way would likely be to perform one step of Gauss-Elimination on A to eliminate the sub-diagonal and then store the inverse of the new main diagonal, the new modified super-diagonal and the constants to preform forward substitution on b. This way each rhs will require one forward substitution to modify b followed by one back-solve. The ...


1

Reiterating some of points made by Tobias and Wolfgang Bangerth: I see no delete to free memory! I highly suggest running valgrind on your code to find any memory leaks. Storing your matrix with two dimensional arrays probably isnt the best idea especially for large systems. Instead I suggest looking into other sparse matrix storage formats, e.g. Compressed ...


1

It is generally true for the $LU$ factorization of any matrix that the bandwidth of the factors equals the bandwidth of the original matrix. In fact, the result is stronger: that the skyline of the factors is contained in the skyline of the original matrix. This is actually not all that difficult to see if you consider the steps one has to take to compute ...


1

If this is a large band matrix, I assume it is sparse and you could well use the sparse solvers in MKL (such as Pardiso or iterative such as conjugate gradients) http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-78889273-7E77-426A-9B5E-23A7C2378D78.htm


1

This does not directly answer your question, but maybe you want to use something like ARPACK, for example scipy.sparse.linalg.eigsh in python.


1

Great observation to see that $u_0$ can be eliminated. Step back and think about the problem for a second. Specifying a Laplace equation fundamentally states that each point is the average of its neighbors. This is commonly visualized as a rubber sheet, and helps me to think about these things. (Poisson is similar w/ more or less stretchy points) When ...


1

While all the present answers are valid solutions to the practical problem, technically the answer to the question in your title (how to reorder variable to minimize bandwidth) is "it's an NP-complete problem". Article: https://link.springer.com/article/10.1007/BF02280884 .


Only top voted, non community-wiki answers of a minimum length are eligible