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.
1,124
questions
5
votes
1
answer
459
views
Converting distance matrix back into original data
Suppose that we have $N$ points, and a distance matrix $D \in \mathbb{R}^{N \times N}$ describing the Euclidean distance among those points. For now, assume that we do not necessarily know how many ...
- 165
6
votes
1
answer
188
views
Fast way of computing the action of a matrix power on a vector
For integer $k>0$, it is well-known that one can use binary exponentiation to evaluate the matrix power $\mathbf A^k$, where $\mathbf A$ is an $n\times n$ matrix.
However, it is not clear to me if ...
- 61
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
9
votes
2
answers
431
views
Efficiently estimating trace of a product of matrices
I have $d\times d$ real-valued matrices $A_1,\ldots,A_k$, $1000<d<4000$, $k\approx 50$, and need to estimate the trace of the following matrix product
$$t=\text{tr}(A_1 A_2\cdots A_k A_k^T \...
- 2,273
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
2
votes
2
answers
112
views
How to find fundamental matrix based on other fundamental matrix and camera movement?
I am trying to speed up some multi-camera system that relies on calculation of fundamental matrices between each camera pair.
Please notice the following is pseudocode. ...
- 121
0
votes
1
answer
65
views
why is $Hx= x-2u^{H}xu$ ? why not $Hx= x(x-uu^H) + 2((u^Hx)u)^2? $
I have some confusion in this diagram
My confusion : why is $Hx= x-2u^{H}xu$ ? why not $Hx=x(x-uu^H) + 2((u^Hx)u)^2?
$
My thinking is that by using pythogoras theorem blue line(vector) denotes ...
- 101
0
votes
2
answers
843
views
Incomplete Cholesky factorization algorithm
I want to implement incomplete Cholesky factorization to precondition, the algorithm I refer from incomplete Cholesky factorization,
...
- 11
1
vote
0
answers
291
views
Trouble inverting complex matrix with numpy and scipy
I have some matrix-valued, complex data $Z(f)$ with $f\in\{f_0,f_1,\dots\}$ and $Z(f_i)$ being a 3x3 matrix. I require the inverse $Z^{-1}(f)$ in my workflow. After encountering some problems with my ...
- 11
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
10
votes
1
answer
500
views
Generalization of eigendecomposition problem
Let $A\in \mathbb{R}^{n\times n}$ and $v \in \mathbb{R}^n$. We recognize $Av=\lambda v$ for some scalar $\lambda$ as an eigendecomposition problem. Suppose $\mu \in \mathbb{R}^n$, and let $\odot$ ...
- 173
4
votes
1
answer
297
views
Is there any way/any python function to calculate the condition number of the roots of a polynomial directly?
I know that NumPy has linalg.cond(A) to find the condition number of a matrix A. But, if I want to find the condition numbers of the roots of a large polynomial ...
- 41
3
votes
1
answer
238
views
Can someone explain the equivalence between Oja's rule and PCA in a simple way?
I have to give a presentation on unsupervised learning in 2 days, and I have to explain/show the equivalence between Hebb's learning rule (or Oja's rule to be more specific) and PCA. The thing is that ...
2
votes
0
answers
110
views
Regularisation of ill-conditioned matrix-vector problem
I have a linear* problem which arises from an integro-differential system, and writes:
$$
(\mathbf{I}+\lambda \mathbf{A})x = b
$$
where $\mathbf{A}$ is a real full matrix, size $n\times n$, but is not ...
- 362
6
votes
1
answer
424
views
PETSc-like library for Julia
I want to build an application for Material Point Method (and probably other meshfree methods too) in Julia and I am looking for library for direct and iterative solvers that can help me with it. One ...
- 163
5
votes
2
answers
584
views
Computational method to compute both the (log) determinant and inverse of a matrix
Suppose I have a square matrix $\mathbf{A} \in \mathbb{R}^{n\times n}$ and a vector $\mathbf{b}\in\mathbb{R}^n$. In my application I need to accomplish two things.
I need to find the solution of the ...
- 183
4
votes
1
answer
134
views
Solving for a single element of a solution of a linear system
I wish to solve a linear system $A x =b$ in which $A$ is dense but not too large, say no larger than $10\times10$. However, I am not interested in the full solution vector $x = [x_0, x_1, \dots]$, ...
- 725
2
votes
0
answers
80
views
Multigrid method: linear solver and modified residual
I am trying to better understand the FAS multigrid algorithm for Euler equation in FV discretization. The usage of the modified residual (the residual with forcing) inside the different cases:
...
- 410
3
votes
1
answer
110
views
Find $x$ that satisfy $(I-A^*A)+x(\frac{A+A^*}{2})\prec0$ using LMI or SDP on Matlab
Given $A\in\mathbb{C}^{n\times n}$, I want to use LMI or SDP to find feasibility of $x>0$ in the following inequality:
$$(I-A^*A)+x(\frac{A+A^*}{2})\prec0,$$
where $D\prec0$ means that $D$ is ...
- 183
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
0
votes
2
answers
133
views
Implementation of $[X, \cdot]$ as an $n^2 \times n^2$ matrix, where $X$ is an $n \times n$ matrix
Let $M_n(\mathbb{R})$ denote the set of $n\times n$ matrices with real entries. I have an $n\times n$ matrix $X\in M_n(\mathbb{R})$, and I would like to implement the linear operator $[X, \cdot] : M_n(...
- 185
5
votes
2
answers
173
views
Optimizing a quadratic form integral over unit sphere
I have an optimization problem, which is to maximize the following integral over the unit sphere:
$$
\max_B \int d\Omega \mathbf{f}^{\dagger}(\theta,\phi) (B^{\dagger} + B) \mathbf{f}(\theta,\phi)
$$
...
- 51
0
votes
0
answers
65
views
Comparing minimas of two different functions
The goal is to find vectors $x_u$ and $y_i$, both of the same length $f=64$, and to do this the following loss function is minimized:
$$\sum_{u, i} (1 + \alpha \cdot r_{ui})(p_{ui} - x_{u}^{T}y_i)^2$$
...
- 41
6
votes
1
answer
249
views
Algorithm for solving systems which are nearly symmetric/adjoint?
I am familiar with Cholesky decomposition and LU factorization for solving systems of linear equations.
I have a problem where I have large sparse matrices (say, 1000x1000 or larger) where only one or ...
- 173
19
votes
3
answers
3k
views
Why do we usually not want the eigenvalues of non-symmetric matrices?
I came across this line in a class note I am reading where it discusses finding eigenvalues of matrices.
In reality we don't go all the way with Arnoldi. We stop at a decent value of 𝑘. Then the 𝑘 ...
- 383
5
votes
1
answer
182
views
Solving absolute value systems
Let $z, b \in \mathbb R^n$, $A \in M_n (\mathbb R)$ and $|z| := (|z_1|, \dots, |z_n|)$. I am searching for an efficient algorithm to solve the absolute value system:
\begin{equation}
z - A |z| = b.
\...
- 111
1
vote
1
answer
690
views
Efficiency of scipy.sparse.linalg.expm_multiply with sparse vs unsparse vectors
From the package scipy.sparse.linalg in Python, calling expm_multiply(X, v) allows you to compute the vector ...
- 185
0
votes
0
answers
86
views
PLASMA usage (Linear algebra routines that supports multithreading)
I have been looking for linear algebra libraries that support multithreading. I have found PLASMA which looks promising. It is from the same group that developed LAPACK.
http://icl.cs.utk.edu/...
- 45
0
votes
0
answers
163
views
Do the projection (with Jacobian) and marginalisation (inversion of matrix and remove a row/column and reinversion) commute?
I try to check the equality or the inequality between 2 Fisher matrices.
The goal is too see if the projection (with Jacobian) and marginalisation (inversion of matrix and remove a row/column and ...
user14831
3
votes
1
answer
394
views
How can I extract the banded or block diagonal part of a sparse matrix in MATLAB?
Given a large sparse (square) matrix in MATLAB, how can I extract the banded or the block-diagonal parts (of fixed size) of it efficiently?
These are useful operations when prototyping and testing ...
- 2,411
0
votes
1
answer
108
views
Efficient solution to linear system involving Kronecker sum in MATLAB
High dimensional finite difference problems often lead to linear systems of the form
$$
A x = b, \qquad A = B_1 \oplus B_2 \oplus \cdots \oplus B_d,
$$
where $\oplus$ denotes the Kronecker sum. $B_i \...
- 1,074
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
1
vote
0
answers
56
views
Number of words that a processor can handle for PMMM
The PMMM which stands for parallel matrix-matrix multiplication essentially accelerates the algorithm of the matrix-matrix multiplication of two matrices $A$ and $B$ both of size $n$ so that $C:= AB$. ...
- 169
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
1
vote
0
answers
83
views
Solution of an underdetermined system stemming from a PDE with Neumann BC
Consider the Poisson's equation in 1D with homogeneous b.c.'s $\mathrm{d} \phi/\mathrm{d} x=0$ with the seven point Laplacian (1 -54 783 -1460 783 -54 1 / 576 on a uniform grid). The resulting system ...
2
votes
4
answers
854
views
Eigenvalue decomposition for a very huge matrix of medical images (such as the pixel physical coordinates of CT images)
I am trying to do eigenvalue decomposition for a huge matrix larger than 788000×788000 for medical image analysis. The matrix is not sparse and every element in the matrix has a real value. And, for ...
4
votes
1
answer
131
views
Solving geodesics on triangular meshes gives negative distances
I have implemented the heat method for geodesics:
https://www.cs.cmu.edu/~kmcrane/Projects/HeatMethod/paperCACM.pdf
When I run it I am getting a solution that, visually, seems correct:
In this image, ...
- 263
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
3
votes
1
answer
324
views
Lanczos algorithm for finding top eigenvalues of a matrix sum
I am trying to find the top k leading eigenvalues of a NumPy matrix (using python dot product notation) L@L + a*[email protected], where $L$ ...
- 203
0
votes
0
answers
486
views
Numerical Range of a matrix in Python
In the mathematical field of linear algebra and convex analysis, the numerical range or field of values of a complex
$n\times n$ matrix $A$ is the set
$$W(A)=\left\{{\frac {{\mathbf {x}}^{*}A{\...
user38211
2
votes
3
answers
194
views
Linear solver recommendation(s) for small problems
I am interested in solving many linear systems $Ax = b$, where $A$ is symmetric positive definite and small (i.e. less than 25,000 rows) --- $b$ will be changing.
We can assume that $A$ arises from ...
- 617
2
votes
0
answers
110
views
Software for solving large systems of linear equations over gf(2)
What available solvers are there for linear equation solver over GF(2) (Boolean), capable of dealing with large sparse systems (in the 10k - 100k variables range)?
- 21
-1
votes
1
answer
57
views
Interpretation of error between Hessian approximation and real Hessian - Quasi-Newton Method
$$ ||I- H_{k}^{BFGS}\nabla^{2}f(x_{k})||_{2}$$
, where $H_{k}$ is the inverse of hessian approximation at each iteration.
I am given this expression to assess the error in Hessian approximation in ...
- 1
-2
votes
1
answer
92
views
Forward Euler Adaptive Step Size Stability
Given
with a generalization using adaptive times-stepping as
then is it still reasonable to assume that to ensure stability of the Euler’s forward method we need the growth factor for all n to be ...
- 1
8
votes
3
answers
1k
views
Accurate Way to Calculate Matrix Powers and Matrix Exponential for Sparse Positive Semidefinite Matrices
I do need to numerically calculate the following forms for any $x\in\mathbb{R}^n$, possibly in python:
$x^T M^k x$, where $M\in\mathbb{R^{n\times n}}$ is a PSD sparse matrix, $n$ can be quite large ...
- 183
2
votes
0
answers
67
views
How to solve this boundary value problem which has more unknown than equation on MATLAB
I need your helps about solving the problem below with MATLAB. I am trying to solve 2D Stress Wave Propagation problem by using FDTD(Finite difference time domain) method on the cylindrical coord. I ...
- 21
1
vote
2
answers
844
views
Diagonalization using LAPACK
Say, we have a Hamiltonian which for simplicity does not mix particle hole sectors. It is just a simple Hamiltonian in real space as shown,
$H=\sum_{ij,\sigma} A(i,j)(c_{i\sigma}^{\dagger}c_{j\sigma} +...
- 13
0
votes
2
answers
408
views
Given a symmetric matrix, is it ok to apply Cholesky decomposition to see if it has negative eigenvalues?
I intend to check the diagonal of L, where A = L'L, for negative elements. However, I don't know if Cholesky is meaningful in theoretical / computational sense if there are some negative eigenvalues.
- 103
0
votes
1
answer
143
views
Solve for large array of PD matrices
I have N matrices that are positive definite, and I have to solve for a M vectors.
As M is large in my case, doing all solves simultaneously using np.linalg.solve ...
- 109
1
vote
1
answer
152
views
Algebraic multigrid as solver and as preconditioner
My question is around the efficiency of AMG. In which case AMG can perform better,as solver or as a preconditioner(for example a Krylov space method as CG)? Assume the case of elliptic pdes.
- 481