Questions tagged [linear-solver]
Referring to methods for solving linear systems of equations.
300
questions
2
votes
1answer
116 views
How to solve for f(A)x=b without GMRES?
How to solve for $f(A)x=b$?
For GMRES, an answer is given in this book chapter: http://link.springer.com/chapter/10.1007%2F978-3-642-58333-9_2.
Ungated version:
https://www.researchgate.net/profile/...
1
vote
0answers
48 views
Residual of Poisson equation with periodic boundaries
I am trying to write a multigrid solver for Poisson's equation, $-\Delta u=f$, on the unit square, $\Omega=(0,1)^2$ with periodic boundaries. My primary source has been Multigrid by Trottenberg, ...
2
votes
0answers
58 views
Solving a huge least squares system of equations when I can only evaluate Ax
I have a situation where I can generate a system of $M$ linear equations for $N$ variables ($N \ne M$). Implicitly this is of the form $Ax=b$ with $A \in \mathbb{R}^{M \times N}$, although I never ...
1
vote
2answers
84 views
How to understand the storage of the Hessenberg matrix of Krylov subspace matrix?
For the Krylov subspace method to solve the large sparse linear system, we first need to generate a subspace Km = span{v,Av,...A^{m-1}v}, which indeed a process ...
0
votes
0answers
32 views
How to implement the gmres method using Householder transformation instead of the Gram-Schmidt?
For Generalized Minimal Residual method GMRES, we usually use the Modified Gram-Schmidt MGS to generate an orthonormal basis of ...
3
votes
0answers
51 views
Levinson Recursion for Non Square Toeplitz Matrices
Given a rectangular Toeplitz Matrix $ H $, how could one solve:
$$ y = H x $$
For instance, $ H $ can be Linear Convolution Matrix of the filter $ h $:
$$ H = \begin{bmatrix}
{h}_{1} & 0 & ...
2
votes
1answer
80 views
Testing a block tridiagonal system of equations
In 1D problems, tridiagonal systems of equations are obtained when we use finite-difference or finite-volumes in a structured mesh. A wide solver is the TDMA algorithm here. In two-dimensional ...
4
votes
0answers
41 views
Nonlinear least squares optimized Jacobian calculation
I have a nonlinear least squares problem, in which I am trying to minimize residuals which can be divided into four classes:
$$
\min_x ||\epsilon(x)||^2 + ||\xi(x)||^2 + ||\delta(x)||^2 + ||s(x)||^2
$$...
4
votes
3answers
148 views
Solving for a vector in a linear system that is both left and right multiplied
I have a linear system where I am given 2 matrices, $A$ and $B$, and 2 vectors, $v$ and $c$, and I need to solve for the vector $x$. $A$ is $n\times n$, $B$ is $n \times n \times n$, and the vectors $...
2
votes
1answer
52 views
What is the standard, extrapolation, and modified version of Richardson iteration method?
I have been studying the iterative methods recently. For classical iterative methods solving $Ax=b$, I have seen that the most simplest iteration method is the so-called "Richardson iteration". But I ...
2
votes
1answer
67 views
What method to solve a sparse complex symmetric (non-Hermitian) system?
I have a sparse system (about 78% of zero entries) that is complex and symmetric (but not Hermitian). The following figure shows the structure of the problem. The off-diagonal blocks are incidence ...
2
votes
1answer
114 views
Why Krylov subspace iterative methods are faster than classical iteration?
This semester, I have been studying the most popular iterative methods, i.e., Krylov subspace iteration methods.
For a large sparse system linear
$$
Ax=b,
$$
where $A$ is nonsingular, I know that ...
1
vote
0answers
29 views
What kind of problem or matrices are suitable for multigrid method?
For Poisson or Convection-diffusion equation as follows:
$$
-\Delta u=f,\qquad u|_\Omega = g.
$$
or
$$
-\Delta u +\vec{w}.\nabla u=f,\qquad u|_\Omega = g.
$$
using FDM or FEM discretization, we can ...
0
votes
0answers
37 views
Why does the initial guess for linear system usually choose by zero vector?
For solving linear system
$$
Ax=b,
$$
using iterative mehods, we often use the terminate criterion as follows:
$$
\frac{\|r_k\|}{\|r_0\|}=\frac{\|b-Ax_k\|}{\|b-Ax_0\|}<eps.
$$where $x_0$ is the ...
2
votes
1answer
50 views
Partially Banded Matrix
I have a somewhat peculiar Jx=R system that I need to solve. The matrix J is 2N -by- 2N. The first N rows have all entries filled. The next N rows are banded in two places, i.e. for the (N+k)th row, ...
0
votes
1answer
59 views
Method to calculate solution of a linear equation system?
I am searching a solution method for the following equation system of equation systems:
Let $A, B \in \mathbb{R}^{n \times n}$ be s.p.d. Matrices and $O$ be the zero matrix of the same size. Further ...
5
votes
0answers
144 views
Symmetric sparse direct solvers in scipy
scipy.linalg.solve, in its newer versions, has a parameter assume_a that can be used to specify that the matrix $A$ is symmetric ...
1
vote
0answers
50 views
Hybrid Ellpack-Itpack (ELL) + COO Sparse Matrix Representation decomposition threshold
Hybrid ELL-COO sparse matrix representation can be done as in the picture, I was looking intensively, however I couldn't find out what is the threshold of decomposing the original matrix into ELL part ...
1
vote
0answers
57 views
implementation for coppersmith matrix multiplication
Is there any online implementation for the coppersmith matrix multiplication I have searched alot but can not find any? and if there is not any why is that Isn't this algotithm much faster than ...
0
votes
0answers
136 views
What algorithm do BLAS and ATLAS use for matrix multiplication?
I have searched and what I understood was that they use the naive one with several memory and cache optimizations.
However, I wanted to know whether they are using the Strassen or the Coppersmith-...
1
vote
1answer
95 views
Lapack symmetric update $B^{-1}AB^{-T}$
Does Lapack have a routine that, given symmetric $A=A^T$ and $B$, computes the symmetric matrix $B^{-1}AB^{-T}$ (while preserving symmetry exactly)?
It would be enough to have this routine for ...
0
votes
1answer
62 views
The proper way to assess the error of Jacobi iteration (for 2D Poisson equation)?
Motivation: I'm using 2D regular grid (it's actually a quadtree but I can still treat it as a finite difference thing if I weight-average the solution over smaller scale cells for the purpose of ...
2
votes
2answers
205 views
Why iterative method: AMG preconditioned PCG is slower than Matlab direct method 'A\b'?
Recently, I have met a question that
a saying goes that for large linear system: iterative methods are required because of memory problem of direct methods.
But when I implement some experiments ...
3
votes
1answer
70 views
Why are the round-off errors when solving the linear system $Ax = b$ of order $\varepsilon_\text{mach} x_j$?
I was reading a paper on arXiv where, in Section 2.4, the authors are discussing the error that arises in the solution of a linear system
$$Ax = b,$$
or, to match up better with the paper,
$$\Phi \...
1
vote
1answer
91 views
Sparse matrix inversion
I have the impedance matrix $Y$, formulated from an electrical network by augmented nodal analysis. The matrix $Y$ is shown as an image to illustrate its feature visually, where all the white blocks ...
0
votes
1answer
137 views
How to set an initial guess for the iterative solver in Comsol?
How to set the initial guess for the iterative solver GMRES or FGMRES for linear problems (Helmholtz equation of RF module) in Comsol?
5
votes
1answer
143 views
Complex differentiation of linear solvers
I have a linear system $$Ax=b$$ which I'm solving approximately, and I need to take the frechet derivative of x with respect to z. Were I solving the problem exactly (either analytically or to machine ...
0
votes
0answers
74 views
Algebraic multigrid for coupled equations
As far as I understand is algebraic multigrid(AMG) a method that was intentionally developed to solve linear systems where every grid point or node has a single DOF. When AMG should now be used for ...
7
votes
0answers
77 views
How to construct an effective preconditioner for this particular problem
A quick introduction to my problem
I am currently developing a method for simulation of water waves in three dimensions based on potential flow theory. The computational bottleneck of the method is ...
1
vote
0answers
45 views
Stability of SVD, Eigendecompositions, and pseudoinverse procedures in modern LAPACK routines
I have proposed an optimisation algorithm which I claim has improved upon the previous algorithm in a number of ways.
One of these claims is that my proposed solution requires no explicit SVD and ...
2
votes
1answer
343 views
A fast way to check if a Matrix is ill-conditioned, and turning it into well-conditioned
I'm running a simulation, and some linear solvers are returning a message of ill-conditioned matrix.
Hence, I'm looking for a fast, easy to implement, method to detect if a matrix is ill-conditioned, ...
6
votes
2answers
116 views
Efficient approach for solving matrix plus diagonal matrix system that varies in time
When solving a system of ODEs, as part of a preconditioner, I get the system
$(A + D(t))x = b(t)$
where $A$ is a sparse matrix and $D(t)$ is diagonal. I'm currently solving this by taking the LU-...
2
votes
1answer
82 views
Assessing numerical error in solving a least squares problem
I have a linear system of the type
$$Ax = b$$
I want to minimise $|b - Ax|^2$.
I know there are different approaches to directly solve the system (Normal equation + Cholesky, QR decomposition, SVD ...
2
votes
0answers
100 views
Richardson's Iteration, Gradient Method and Spectral Radius
Richardson's iteration introduce a scalar $\alpha$ to the update formula:
$$ \textbf{x}^{(k+1)} = \textbf{x}^{(k)} + \alpha \textbf{r}^{(k)} $$
And compute $\alpha$ by minimizing the spectral radius:...
2
votes
1answer
64 views
Right-preconditioning and fixed point linear iterations
Given a linear system $A\textbf{x}=\textbf{b}$, we can express it into the easier-to-solve right-preconditioned form:
$$ AM^{-1}\textbf{y}=\textbf{b}, \quad \textbf{y}= M\textbf{x} $$
On the other ...
9
votes
4answers
434 views
fastest linear system solve for small square matrices (10x10)
I am very interested in optimizing the hell out of linear system solving for small matrices (10x10), sometimes called tiny matrices. Is there a ready solution for this? The matrix can be assumed ...
3
votes
0answers
65 views
Numerical analysis, pivoting and incomplete LU decomposition
When doing LU decomposition, the algorithm will break down if any of the diagonal element $x_{ii}$ is zero. Therefore, we can use pivoting on the matrix such that $x_{ii}$ is no longer zero. That is ...
6
votes
1answer
108 views
numerical solution of an under-determined linear equation in high dimensions
I need to solve a linear regression problem $$Ax=y$$ which is hugely underdetermined. I have around $10^6$ features but only $10^3$ equations. So $A$ is a $1,000\times 1,000,000$ matrix and $y$ a ...
2
votes
2answers
163 views
Asymptotic Complexity of Gaussian Elimination using Complete Pivoting
I would like to know the algorithm asymptotic complexity with Complete Pivoting. With partial pivoting, it is known to be $O(n^3)$.
Is it the same for complete pivoting?
1
vote
1answer
94 views
Conjugate gradient - ill-conditioning and numerical tolerance
I would like to solve system $Ax=b$, where $A$ is SPD, but very ill-conditioned ($\text{cond}(A)>10^{11}$). I am interested in using UNpreconditioned version of the conjugate gradient method.
Is ...
4
votes
1answer
2k views
Which C++ linear algebra library is probably the fastest on solving huge sparse [square matrix] linear system?
I am developing a 2D CFD solver for fluid-particle interaction. To solve Navier-Stokes equations on a grid of size $10000\times 10000$ cells (or >1 million cells), a large linear system $Ax=b$ with $A$...
2
votes
1answer
341 views
Memory and time requirements of the scipy sparse spsolve
I have a system of fairly large set of linear equations (approximately 30K equations). I am using scipy.sparse.spsolve to solve these equations. Initially, I tried ...
0
votes
1answer
36 views
Can LINCS algorithm be used for colliding molecules?
Supposing that one molecule is static and one is dynamic,
can the dynamic one be solved with LINCS for its shape (angle, bond length) constraints and also keep collisions with static molecule off, ...
2
votes
1answer
96 views
Is steady linear elasticity inherently ill-conditioned?
Compared to the transient PDE for linear elasticity, the steady equations appear to less well-conditioned. Are they inherently ill-conditioned without the transient term?
The condition number for ...
2
votes
1answer
127 views
Condition number of matrix and effects of round off errors
In my numerical linear algebra class, I learned that for some matrices, it could have an element that is a very small number that is approximately 0 (and many orders of magnitude different from all ...
3
votes
1answer
154 views
Derivatives of Approximate Matrix inverses
I am cross posting this question to the mathermatics stack exchange. please find it either at this link, https://math.stackexchange.com/q/2952989/430980, or below:
I have a question concerning the ...
8
votes
0answers
182 views
Imbalance of variables in Mixing Newton's method and Linear solver for a Non-linear system
Problem
Solving a non-linear system of equations.
The number of variables is the same as the number of equations.
When I fix a set of variables (say $\vec{y}$) and keep another set free (say $\vec{...
2
votes
2answers
128 views
Memory/speed tradeoff for many small matrix inverses
Problem
In the case of a finite element code, I have many small (order of 30x30) matrix inverses (or LU factorizations), one per finite element. These matrix inverses never change and must be applied ...
0
votes
1answer
63 views
How to access solution to linear system in PETSc?
I have just started with PETSc hence it might seem like a very stupid question but I couldn't find any answer in manual.
After Calling KSPSolve, where can I access the soluition for my linear system?
...
3
votes
1answer
151 views
How many operations are needed for LAPACK's zgesv to solve a linear system?
I have a linear system of complex numbers. I am using LAPACK' zgesv (actually I am using intel MKL LAPACKE, but I am assuming the algorithm is the same). No assumption can be made about the system.
I ...
