Questions tagged [linear-algebra]
Questions on the algorithmic/computational aspects of linear algebra, including the solution of linear systems, least squares problems, eigenproblems, and other such matters.
213
questions with no upvoted or accepted answers
13
votes
0
answers
707
views
Fast Eigenvalue and SVD Solver for Structured Matrices
I am looking for a fast Eigenvalue and SVD solver for small dense structured matrices (Hankel and Toeplitz). I have searched for efficient implementations in libraries like MKL but I am not able to ...
- 241
9
votes
0
answers
229
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{...
- 163
9
votes
0
answers
293
views
Wanted: sequences of linear systems for recycling Krylov solver analysis
In the solution of sequences of linear systems
$$A_ix_i=b_i\quad\text{for}\quad i=1,2,\dots$$
with Krylov subspace methods, data can be recycled from already solved linear systems in order to speed up ...
- 231
8
votes
0
answers
281
views
Stable alternatives to "condition number"?
A number of numerical problems are easy to solve when condition number $\kappa$ of the problem is low. For instance, conjugate gradient descent complexity scales as $O(\sqrt{\kappa})$.
However
"...
- 2,273
8
votes
0
answers
426
views
Eigenvalue with largest imaginary part
Iterative eigensolvers such as ARPACK, give the option to find a subset of the eigenvalues which have the largest imaginary part. My question is how do these algorithms work.
As I understand it, ...
- 243
8
votes
0
answers
769
views
What is the source of the error in the Sherman-Morrison formula application?
The Sherman-Morrison formula
$$ (A+uv^T)^{-1} = A^{-1} - \frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u} $$
results in small errors in relation to the standard matrix inverse operation after each application, ...
- 180
7
votes
0
answers
84
views
Choice between using Moore-Penrose inverse and G2 inverse
Moore-Penrose inverse for an arbitrary matrix $X\in \mathbb{R}^{n \times p}$ is defined by a matrix $X^\dagger$ satisfying all of the Moore-Penrose conditions, namely
\begin{align}
(1) \;\;\;& XX^\...
- 173
7
votes
0
answers
135
views
Can we sparse solve a few eigenvalues specified by index range?
I need to solve a few eigenvalues of a large sparse matrix specified by their index range. These indices are according to the whole eigenspectrum sorted in algebraic (not absolute value) ascending ...
- 303
7
votes
0
answers
333
views
Implementation of Lanczos method that returns tridiagonal matrix
The Lanczos method can be used to obtain extremal eigenpairs of sparse symmetric or hermitian matrices. I know there are several implementations of the Lanczos method (as well as Arnoldi, Davidson, ...
- 171
7
votes
0
answers
326
views
Compute sparsity pattern of $A^2$
Suppose we have a sparse matrix $A$. Is there any way to compute just the sparsity pattern of $A^2 = A \cdot A$ (I do not actually need to know what exactly the nonzero value are) faster than to ...
- 71
6
votes
0
answers
508
views
Is there any catch on using `zgemm3m` vs regular `zgemm`?
I've just (to my embarrassment) encountered a BLAS-like extension of a matrix-matrix product subroutine gemm in Intel MKL: gemm3m...
- 8,602
6
votes
0
answers
134
views
What strategies / decompositions would be useful to solve the following linear system repeatedly if I only care about time to solution?
Let $A\in \mathbb{R}^{n\times n}$ be symmetric positive semidefinite, and $B\in \mathbb{R}^{n\times n}$ be symmetric positive definite. Suppose $B$ is block diagonal so it is easy to invert. (We ...
- 1,000
6
votes
0
answers
140
views
How does an unpivoted QR fail to reveal rank?
An unpivoted QR factorization produces a triangular factor $R$. A rank-revealing QR factorization is typically done with column pivoting. My question is, how does an unpivoted QR factorization fail to ...
- 4,480
6
votes
0
answers
193
views
Find the solution of linear equation using Wiedemann/ Krylov method
Let given
$M =$
1 0 1
0 1 1
1 1 1
and
$b =$
1
0
1
How to find the solution $x_3$ where $x=${$...
user8264
5
votes
0
answers
135
views
Why are fast Givens rotations mentioned so little in the recent literature?
Fast Givens rotations seem like a nice upgrade from standard Givens rotations. They take fewer multiplications to apply and avoid the calculation of a square root.
They are, however, given very little ...
- 181
5
votes
0
answers
201
views
Inverse problem with uncertain forward operator
Suppose I want to solve a linear inverse problem. In this example we take a convolution with the kernel:
$$\frac{1}{(y^2+z^2)^{3/2}}$$
We only take a fixed $z$ for the computation and convolve with ...
- 725
5
votes
0
answers
144
views
Explanation of subspace strategy regarding CG described in Golub's book
I was wondering about the last paragraph in Matrix Computations (4th edition) by Golub, Chapter 11 (11.3.3), specifically his explanation of subspace strategy for Conjugate Gradient.
Note that in ...
- 161
5
votes
0
answers
97
views
Trying to understand splitting-based iterative method for 2D Laplace problem
I am trying to understand the theory behind a splitting based iterative method which uses the incomplete Cholesky factorization. Before giving the specific details, let me first give the problem ...
- 51
5
votes
0
answers
130
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 & ...
- 332
5
votes
0
answers
87
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
$$...
- 1,048
5
votes
0
answers
90
views
Improving convergence of Jacobi iteration to Schur form
I'm using SIMD processor arrays to compute the eigen-decomposition for large numbers of small (up to $32\times 32$) matrices. For assorted technical reasons, Jacobi iteration maps well to the SIMD ...
- 51
5
votes
0
answers
158
views
Multi-matrix orthogonal basis problem
Suppose we are given a set of symmetric, positive definite matrices $A_1,A_2,\ldots,A_k\in\mathbb{R}^{n\times n}$. Is there any numerical method or reduction to a known problem (e.g. eigenvalue ...
- 1,341
5
votes
0
answers
240
views
Preconditioning technique for large sparse non-hermitian matrix
I am attempting to solve a computational acoustics problem that involves solving an underlying sparse matrix. The size of the problem varies with grid size (3D) and fill-in's obviously make direct ...
5
votes
0
answers
775
views
Perron-Frobenius theorem on general real symmetric matrices
From the Perron-Frobenius theorem, it might be concluded that the spectral radius is the largest eigenvalue for positive matrices, ie, matrices with strictly positive entries. In other words, the ...
- 1,663
4
votes
0
answers
83
views
Can you left-precondition least squares?
Suppose I want to solve an overdetermined linear least squares problem
$$
x = \operatorname*{argmin}_{x\in\mathbb{R}^n} \| Ax - b\|^2
$$
where $A \in \mathbb{R}^{m\times n}$ has full column rank. I ...
- 350
4
votes
0
answers
100
views
Decomposing a banded matrix
Suppose we have a linear algebra problem with a banded matrix A which has nonzero entries on the main diagonal, two nearest sub-diagonals, and two other sub-diagonals (such band structure often arises ...
- 2,440
4
votes
0
answers
182
views
Stable iterative solver for complex symmetric linear systems
I am interested in the iterative solution (preferably Krylov-type solvers) of a problem $\boldsymbol{A}x=b$, with $x,b\in\mathbb{C}^{n\times1}$ and $\boldsymbol{A}\in\mathbb{C}^{n\times n}$. $\...
- 141
4
votes
0
answers
91
views
Unstable Algorithms which become stable when hardware provides Kulisch exact dot product instruction
In John Gustaffson's book The End of Error, he discusses Ulrich Kulisch's exact dot product, which (in double precision) requires a 2100 bit fixed point register which rounds only once after the ...
- 2,135
4
votes
0
answers
99
views
Check if LinearOperator is symmetric
I have a scipy.sparse.linalg.LinearOperator object. I'd like to check if its associated matrix is symmetric without actually instantiating the matrix in the most ...
- 203
4
votes
0
answers
74
views
Efficient computation of marginalized multivariate normal likelihood
In general,if we know that the marginal Gaussian distribution for some variable $\textbf{x}$ and a conditional Gaussian distribution for some $\textbf{y}|\textbf{x}$ of the forms:
$$p(\textbf{x}) = \...
- 41
4
votes
0
answers
76
views
Optimize linear equation using inner products and subject to L1 norm
I have a linear system of the form $A x = b$ where $A$ and $b$ are known, $A$ is "square", and $\lvert b \rvert_1 = \lvert x \rvert_1 = 1$. Unfortunately, I am working in a framework that ...
- 245
4
votes
0
answers
543
views
Fastest matrix library for Android (with GPU is possible)
I was working on an Android app that requires some linear algebra with matrices. The matrices will be somewhat medium-sized as they are not too small or too big. I was originally using jBlas because ...
4
votes
0
answers
100
views
Block matrix and DSYRK
I want to compute the matrix
$$
A = \sum_{i=1}^N v_i v_i^T
$$
where each $v_i$ is a given vector of length $2500$, so that $A$ is $2500 \times 2500$, and my $N$ is about 2 million. Rather than call ...
- 1,048
4
votes
0
answers
135
views
How to stochastically estimate the trace of a matrix?
Specifically, the diagonal elements (can possibly both positive and negative) of the matrix can be computed efficiently but the total number is large ($\mathcal O(10^{18})$).
My first thought about ...
- 61
4
votes
0
answers
1k
views
Fixing a near singular covariance matrix
Given a near singular covariance matrix, the standard method of 'fixing' it seems to be to add a small damping coefficient $c>0$ to the diagonal, which serves to bump all the eigenvalues up by this ...
- 461
4
votes
0
answers
214
views
Galerkin FEM error when using even number of elements
Intro:
I am developing a Galerkin FEM code in matlab - starting small with a simple 1D ODE. The equation I'm trying to solve is: $$a u_x = cos(x), x \in [0, 2\pi]$$ Which has a known exact solution $$...
- 674
4
votes
0
answers
103
views
Tracking for two meshes
I'm dealing with physical simulation (position based dynamic). Now I'm trying to realize tracking for two meshes. In order to explain what does it mean, let's assume that we have two similar meshes: ...
- 289
3
votes
0
answers
112
views
Stochastic power iteration for generalized eigenvalue problems?
Suppose $\mathbf{x}$ is a random variable in $n$ dimensions, and $u$ is a vector. How can I estimate the following quantity in an online fashion?
$$f(x)=\max_{\|u\|=1} \frac{ E\left[\langle u\cdot x\...
- 2,273
3
votes
1
answer
143
views
Estimating the spectral radius when applying the method of lines
Some time integrators, notably the Runge-Kutta-Chebyshev method, implemented in the RKC code from Sommeijer & Verwer, gives the user an option to provide a callback with an estimate of the ...
- 544
3
votes
0
answers
117
views
When would one choose un-pivoted $LDL^T$ instead of $LL^T$ for a Positive Definite Matrix?
Background
I am decomposing (and then operating on) a symmetric positive definite matrix (a covariance matrix) as part of a larger system. Dimensions range from O(10) to O(300). The covariance ...
- 792
3
votes
0
answers
71
views
Approximating eigenvalues of DPR1 matrix with special properties
In my application, I have a sum of diagonal $A$ and rank-1 $B$
$$T=\underbrace{\text{diag}(1-2\alpha h+2\alpha^2 h^2)}_A + \underbrace{\alpha^2 hh^T}_B$$
Where $h$ is a vector $\in \mathbb{R}^d$ with ...
- 2,273
3
votes
0
answers
77
views
Householder Vector algorithm in Golub and Van Loan
(This is repost of a question first asked on Mathematics. Hopefully there are more people here who have a copy of Golub and Van Loan to hand)
In the 4th edition of "Matrix Computations", ...
- 181
3
votes
0
answers
104
views
Compute orthogonal complement using BLAS / LAPACK
Is there a fast method to compute an orthogonal complement of an arbitrary matrix $U\in\mathbb{R}^{m \times n}$ in BLAS / LAPACK?
Specifically, I want any matrix $V\in \mathbb{R}^{m \times (m - \text{...
- 799
3
votes
0
answers
699
views
Jacobian Matrix of 2D element mapped to 3D
Note: I previously posted this question to MathStackExchange, but got no attention there. So I'm rewritting and trying over here.
Problem summary
Given a common¹ set of shape functions defined at ...
- 39
3
votes
0
answers
746
views
Compute Nullspace of Sparse Matrix
I am computing the nullspace of a sparse rectangular $m$ x $n$ matrix $A$, where $m$ << $n$. I do this by computing the QR decomposition of $A^T$ and extract the $n-m$ right-most columns of the ...
- 661
3
votes
0
answers
102
views
Numerical calculation of the Berry connection
I'm doing some numerical calculations involving Hermitian matrices, and derivatives of the eigenvectors.
Essentially, I have an n x n, Hermitian matrix H(x), which is dependent on some continuous ...
3
votes
0
answers
100
views
Is the matrix exponential and the Jordan canonical form actually useful for solving differential equations?
All of my yearlong graduate-level Linear Algebra course notes from my professor—an algebraist/representation theorist—shows his love for the exponential map $e^A$ and the Jordan canonical form—and one ...
- 31
3
votes
0
answers
361
views
What is the fastest algorithm for computing log determinant of a PSD matrix? (All possible PSD matrices)
I am diving into some literature to understand which is the best algorithm for computing the log-determinant of a PSD matrix. More generally, I am interested in a list of resources to read, which ...
- 131
3
votes
0
answers
47
views
Subspaces for Iterative methods
In the original paper of Conjugate Gradients, the authors mention that if we pick the canonical basis $\{e_1,e_2,\ldots,e_n\}$, to obtain A-orthonormal vectors, we end up with the Gaussian elimination ...
- 171
3
votes
0
answers
665
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-...
- 41