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.

Filter by
Sorted by
Tagged with
0 votes
0 answers
88 views

Is there a fast matrix-free inverse power iteration?

Problem: I want to solve the eigenvalue problem $$x=Ax$$ to the eigenvalue $1$ for a large matrix (roughly $N^3\times N^3$ and $N$ ranges from 10 to 100) where $A$ is stochastic (i.e. all entries are ...
41 votes
5 answers
35k views

How does the MATLAB backslash operator solve $Ax=b$ for square matrices?

I was comparing a few of my codes to "stock" MATLAB codes. I am surprised at the results. I ran a sample code (Sparse Matrix) ...
4 votes
1 answer
211 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 ...
10 votes
5 answers
897 views

How to efficiently compute the determinant of a matrix with unknown diagonal entries?

I would like to ask Python to compute the determinant of a large symmetric matrix where all off diagonal entries are known. The diagonal entries could vary. Since I need to compute the determinant ...
3 votes
1 answer
150 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 ...
2 votes
1 answer
70 views

Computing the Fiedler vector of a large, sparse graph

I have a sparse, undirected and unweighted graph $G$ of size $n$, with $n$ on the order of say several million. I would like to compute the Fiedler vector $f$ of $G$, which is the eigenvector ...
0 votes
0 answers
31 views

Flexible Conjugate Residual

If we want to use variable preconditioning in Conjugate Gradient, we can replace the Fletcher–Reeves by the Polak–Ribière formula (https://en.wikipedia.org/wiki/Conjugate_gradient_method#...
4 votes
1 answer
152 views

Numerically stable computation of $x^T A x$

We have a large sparse symmetric positive-definite matrix $A \in \mathbb R^{N \times N}$ and a vector $x \in \mathbb R^N$. How do I practically compute the inner product $x^T A x$ when the matrix $A$ ...
1 vote
1 answer
68 views

Fill-reducing ordering for computing the matrix product $A^T A$?

I have found many libraries for reducing filling when dong Cholesky factorisation on sparse matrices. However, I want to do fill-reduction for a different reason - given a $m\times n$ matrix $A,$ I ...
0 votes
1 answer
67 views

Weird runtime behavior of `scipy.linalg.solve_triangular` and `trtrs`

I want to understand the time complexity of scipy.linalg.solve_triangular, which calls trtrs from LAPACK under the hood, so I ...
1 vote
0 answers
91 views

How to implement boundary conditions for the Thomas algorithm

For my variable $U(t,x)$, I have implemented the thomas algorithm with $U_j^i$: $$ a(x)U_{j-1}^{i+1}+ b(x)U_j^{i+1} + c(x)U_{j+1}^{i+1} = d(x)U_j^{i} $$ Then $\textbf{A}$ is a tridiagonal vector with ...
0 votes
0 answers
67 views

Best transform of matrix to make it efficient for shift-then-invert?

I am using ARPACK to find the smallest eigenvalue of a matrix. I use the shift and invert method. That is, looking for the largest eigenvalue of $$ (A-\sigma I)^{-1}. $$ However, I do not know $\sigma$...
4 votes
0 answers
189 views

Why for $A^T A$, it is faster to computer the eigenvalues of its inverse than itself?

I have written the following code in MATLAB. I also observe the same effect in Arnoldi iteration in ARPACK in C. ...
0 votes
0 answers
92 views

Computing smallest singular value of a matrix with explicit error control?

[Also posted here: https://mathoverflow.net/q/464433/] Many good algorithms are out there computing truncated SVD: https://mathoverflow.net/q/161252. I am trying to implement some codes to find the ...
1 vote
0 answers
71 views

Ways to speed up Lanczos algorithm when we have a very dense cluster of many eigenvalues?

Lanczos algorithm can be used to find the largest/smallest eigenvalues of matrices. I am trying to find a good library in C/C++/Rust for finding the smallest singular value (or eigenvalue). I have ...
1 vote
2 answers
89 views

Efficient projection onto the kernel of a matrix

Suppose I have a positive semidefinite matrix $M = \sum_i^N A_i^T A_i$ where each $A_i$ is a fat matrix of shape (m,n) and $m << n$, we can also assume that $A_i$ is full rank(but the stacked ...
13 votes
8 answers
3k views

Real-world applications of eigendecomposition?

Cross-posted on Math.SE Are there real-world applications that call specifically for eigenvalues rather than singular values? I often see eigendecomposition used as "poor-man's SVD" For ...
0 votes
0 answers
60 views

Help with inferring Network topology from Spectral templates

I am trying to use matlab and YALMIP to solve a graph learning problem of recovering eigenvalues from the eigenvectors of the covariance of sampled graph signal data. This is to implement the ...
5 votes
2 answers
579 views

$\mathbf{UDU}^\top$ decomposition routines in LAPACK/Eigen?

I would like to compute the decomposition of a real symmetric positive definite matrix $\mathbf{A} = \mathbf{UDU}^\top$. LINPACK seems to have it as DSIFA, but I ...
3 votes
2 answers
179 views

what preconditioner for incompressible hyperelasticity in 3d (similar to stokes equation?)?

I am working on modeling incompressible elasticity at finite strains. $$ \mathrm{Div} \boldsymbol P = \boldsymbol 0, \quad \boldsymbol X \in \Omega_0 \subset \mathbb R^3, \\ J = 1, \qquad \boldsymbol ...
0 votes
0 answers
22 views

dual svm square hinge loss

Let $x_1,\dots,x_n\in \mathbb{R}^n$, $y_1,\dots,y_n\in \{-1,1\}$, $\lambda \ge 0$ and $K$ be the invertible Gram matrix $K=(x_i\cdot x_j)_{ij}$. Consider $$ (P) \qquad \qquad \min_{a\in \mathbb{R}^n} \...
0 votes
0 answers
60 views

Singular Matrix Error in Incomplete LU Decomposition

I’m currently working on solving the following PDE: $$\begin{equation} -(\mu_x \frac{\partial^2 u}{\partial x^2} + \mu_y \frac{\partial^2 u}{\partial y^2}) = f(x, y)\end{equation}$$ Where a right hand ...
5 votes
0 answers
116 views

Name for vectors in a Krylov space but not the preceding one

It seems to me that a useful concept to define when studying Krylov subspace methods is the idea of a vector $v$ that belongs to a Krylov subspace $\mathcal{K}_{n+1}(A,b)$ but not to the preceding one ...
0 votes
0 answers
59 views

Solving AU = F using linalg.cg results in 0 iterations

I am working on solving the following PDE: $$\left(\mu_{x}\frac{\partial^{2}u}{\partial x^{2}}+\mu_{y}\frac{\partial^{2}u}{\partial y^{2}}\right)=f(x,y) \tag 1$$ Which is then discretised: $$- \mu_{x} ...
1 vote
1 answer
86 views

accuracy problem for a null space calculation on a sparse rectangular matrix

I have been using the QR-based approach on this link to find the null space of rectangular matrices, and possibly sparse matrices, that emerge as a result of some coupling conditions of different ...
0 votes
1 answer
75 views

Finding the Vector $v$ for a Given Householder Matrix Transformation of Non-Collinear Vectors $a$ and $b$

Consider a vector $v$ in $\mathbb{R}^{n\times1}$. The Householder matrix is defined as follows: $$H(v)=I-\dfrac{2vv^T}{v^Tv}.$$ It can be demonstrated that $H(v)$ is symmetric and orthonormal. The ...
1 vote
1 answer
88 views

Are there any established direct eigensolvers for sparse hermitian matrices?

I have experience with LAPACK (direct solvers) and ARPACK (sparse iterative solvers), but are there any sparse direct solvers? I am concerned more with preserving space than with fast solutions. ...
12 votes
2 answers
2k views

Why are systems with clustered eigenvalues easy to solve?

I came across the following slide by Theo Diamandis & Zachary Frangella on what makes the linear system $Ax=b$ easy to solve using the conjugate gradient method. Transcription: CG converges ...
2 votes
0 answers
42 views

How to use a preconditioner estimated from a subset of data?

Suppose I'm solving $Ax=b$ using row-action method like Kaczmarz for $m\times n$ matrix A with $m\approx \infty$ and have $H_k=\frac{1}{k}A_k^T A_k$ which is an estimate of the Hessian obtained from ...
3 votes
1 answer
187 views

Solving underdetermined Lyapunov equation?

I'm solving the following for $X$ with $A,B$ singular positive semidefinite matrices. $$AX + XA = B$$ Because $A$, $B$ are singular, standard Lyapunov solver fails However, if I heuristically skip ...
1 vote
1 answer
89 views

enough conditions to check that a matrix doesn't have Cholesky factorization while factorizing it

I wrote this code to find Cholesky factorization of a symmetric positive definite matrix in MatLab: ...
0 votes
1 answer
201 views

Problems on the algebraic theory of sparse matrices

I have finished testing basic large densely parallel matrix multiplication on 4 gpu's ,and have done work on TSLU and TSQR on cpu's based on mpi. I am going to continue working on the theory of ...
3 votes
0 answers
114 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\...
0 votes
0 answers
18 views

How do you determine the Mott-Insulator to Superfluid transition in the Bose Hubbard System

I am doing some simulations on various systems expressed in 2nd quantization and one of the points of interest of mine was Phase transitions in the Bose-Hubbard model $$ H = \sum_{k} \{ t_k(b^\dagger_{...
2 votes
1 answer
1k views

Calculate partial trace of an outer product in Python?

I have a python implementation of calculating the partial trace over select dimensions. ...
1 vote
1 answer
67 views

Dyadic product of a Gaussian vector (in Convolution kernel)

I was reading about the 2D Gaussian-blur convolution kernel. A Gaussian vector is a vector where the elements follow a Gaussian distribution. In my case, the kernel is symmetric and hence we take a ...
11 votes
1 answer
3k views

What is the state of the art algorithm for diagonalizing real symmetric matrices?

There are many methods for diagonalizing matrices, probably the most widely used is the combination of Householder transformations and the QR algorithm. Is there any superior method for diagonalizing ...
4 votes
1 answer
331 views

How to efficient solve $e^{-tA} x =b$, where A is a very sparse matrix

I am going to solve an equation containing an exponential matrix $e^{tA} x =b$, which can be obtained naturally through $x=e^{-tA} b$. A is a 1million $\times$ 1 million matrix with stores 7.15 ...
0 votes
1 answer
287 views

Questions on the theory of distributed numerical algebraic computation

I'm trying to build a pure python distributed numerical algebra computation kernel based on GPU. but after I've learnt most of the software engineering, I realise that I'm seriously lacking in ...
0 votes
0 answers
46 views

recommendation on some papers/books about frontal solver used in FEM

I'm reading a program about computational plasticity, this program use frontal solver to solve the program, but I'm not familiar with frontal solver even after reading some papaers, so could you ...
16 votes
3 answers
16k views

What is the fastest way to compute all eigenvalues of a very big and sparse adjacency matrix in python?

I'm trying to figure out if there is a faster way to compute all the eigenvalues and eigenvectors of a very big and sparse adjacency matrix than using scipy.sparse.linalg.eigsh As far as I know, this ...
0 votes
1 answer
93 views

Computing discrete laplacian matrix for mesh fairing

I asked this question on the math stack exchange and got an answer, but I am just as utterly confused as before. My fundamental goal is to actually construct the matrix, that is, a series of steps I ...
4 votes
1 answer
117 views

Numerical continuation of all eigenvalues of small, dense matrices

Consider a one-parameter family of matrices $A(q)\in \mathbb{R}^{n\times n}$, $q\in\mathbb{R}$. For my applications, $n$ is typically between $5$ and $50$ and $A(q)$ is generally dense, so direct ...
0 votes
1 answer
49 views

How do you build a polyharmonic discrete system?

Polyharmonic equations, to my understanding, are defined as: $$\Delta ^k u = 0$$ i.e. one repeatedly applies the laplace operator to the function a certain number of times and the result must be 0. ...
5 votes
2 answers
1k views

CHOLMOD condition number estimate

The CHOLMOD library provides a CHOLMOD_rcond function that estimates the reciprocal condition number (in the one norm) of a symmetric positive definite matrix from ...
38 votes
10 answers
9k views

stupid + stupid = brilliant in scientific computing

I'm interested in examples of very effective methods in scientific computing that are the sum or naive combination of very ineffective or bad ones.
0 votes
0 answers
47 views

Different computational approaches to show that the conjecture is true

It seems that proving the statement given in the question https://math.stackexchange.com/questions/4601532/a-pen-and-paper-proof-for-a-matrix-implication is difficult analytically. I was thinking what ...
4 votes
1 answer
1k views

Why are all eigen solvers iterative?

I have small dense square matrices for which I would like to compute the inverse by singular value decomposition, or equivalently solve the eigenvalue problem. While there are many direct methods for ...
2 votes
1 answer
175 views

Solution of linear system doesn't work, in parallel

I'm solving $Ax = b$ with PETSc, $A$ sparse and asymmetric. I'm using BCGS or FGMRES or TFQMR as a solver, and ILU as a preconditioner. When I use 1 core, everything works as expected. But with 8 ...
0 votes
1 answer
69 views

Solve bivariate polynomial system

Given a bivariate polynomial system with variables $(x, y, z)$ like (1) $ f_1 = x * a_1 + y * a_2 + z * a_3 = 0 $ (2) $ f_2 = x * a_4 + y * a_5 + z * a_6 = 0$ (3) $ f_3 = x^2 + y^2 - 1 = 0$ how do I ...

1
2 3 4 5
23