Stack Exchange Network

Stack Exchange network consists of 175 Q&A communities including Stack Overflow, the largest, most trusted online community for developers to learn, share their knowledge, and build their careers.

Visit Stack Exchange

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.

3
votes
0answers
79 views

Computing Algebraic Riccati inequality

In my Robust model control, I have got a couple of quadratic Riccati inequalities which need to be solved numerically on MATLAB. My question, Is there function on MATLAB can solve quadratic Riccati ...
5
votes
2answers
237 views

Efficiently removing projection to subspace without having an orthogonal basis

I have a number of vectors $v_1, …, v_n$ and another vector $w$, that all are linearly independent, but not orthogonal. Let $V := \mathop{\text{span}}(v_1, …, v_n)$. I need to remove $w$’s projection ...
1
vote
0answers
86 views

QR via Householder: less computationally complex variants?

I'm a probabilist and need to do a few computations for a rather big linear least squares problem, so I'm trying to optimize the computation as far as is feasible to me. In computing the QR ...
0
votes
1answer
107 views

Role of non-hermitian coefficient matrices in the discretization of self-adjoint operators

What is the role of non-symmetric coefficient matrices in the solution of partial differential equations with self-adjoint operators? Particularly, I'm thinking about time-propagation of a linear ...
1
vote
0answers
53 views

Closed-form Jacobian of se3 element w.r.t. 6-dof motion

Let $A$ and $B$ be two rigid transformations in 3D space that transform things from global to local coordinates. Let their relative transformation be expressed by $W=A*B^{-1}$. $W$ can also be ...
1
vote
1answer
104 views

Finding the lowest $n$ eigenvalues of a band-diagonal Matrix

I have a real sparse matrix of the form $$ \left( \begin{array}{ccc} h_{11} & h_{12} & 0 & h_{14} & & & \\ h_{21} & h_{22} & h_{23} & 0 & h_{25} & & ...
1
vote
3answers
126 views

Solving a small non-symmetric, non-diagonally dominant, and non-sparse system

I want to solve a small (20 $\times$ 20 up to 30$\times$30) system which is not symmetric, not diagonally dominant, and not sparse. Each row contains a modified form of the Legendre coefficient of a ...
3
votes
0answers
150 views

Eigenvalues with absolute values close to 1

I have a generalized eigen-problem: $A\psi = \lambda B \psi$ with $A$ and $B$ are large (>1000) complex non-Hermitian matrices. I know that eigenvalues with largest and smallest absolute values can be ...
2
votes
0answers
158 views

Does applying the Newton-Raphson iteration for matrix reciprocal refine a matrix inverse from LU/GE?

This is a follow-up to this answer. Suppose you have a possibly very ill-conditioned matrix $A$, and you compute its inverse with LU/GE to get $X_{\text{lu}}\approx A^{-1}$. The Newton-Raphson ...
1
vote
1answer
160 views

What are the prominent algorithms for solving systems of linear inequalities?

Why, you may ask? The Python sympy symbolic library provides solutions to single, linear Diophantine equations in terms of parametric variables. For instance, ...
5
votes
3answers
803 views

Least Squares: Numerically, is solving normal equations okay for nice matrices?

I have to solve a least squares problem: $$ x=\arg \min\|Ax-b\| $$ where $A$ is a $m\times n$ matrix, $m>n$, and $b\in\mathbb{R}^m$. I always thought that doing this via QR factorization is ...
4
votes
2answers
637 views

Smart way to multiply 3 matrices

I have a quantum mechanics simulation where I need to multiply three matrices that look like this: $$\rho(t_1)=U^\dagger \rho(t_0) \, U$$ where $U^\dagger$ is the hermitian conjugate of $U$. This ...
1
vote
0answers
156 views

How to obtain projections from sinogram in ART reconstruction technique?

I'm kind new in the Computed Tomography field and I'm trying to understand and implement ART technique. Said it, I started to read the book The Mathematics of Medical Imaging - A Beginners Guide by ...
2
votes
2answers
1k views

Fastest way to perform element-wise multiplication on a sparse matrix

I have two large-ish matrices (~100K cols x ~100K rows). They are sparse and symmetrical (about 0.1% of them values are non-zero). I want to do element-wise multiplication between them. Also, I ...
0
votes
1answer
146 views

How to compute matrix representation of $\hat{y}\frac{\partial}{\partial x}$?

I have a 2-dimensional system which I would like to solve numerically (I'm using finite difference method right now), and its an eigenvalue problem. I have a term that looks like $H\psi(x,y) = [-\frac{...
4
votes
1answer
507 views

Thomas algorithm for 3D finite difference

For 1D finite difference, the resulting linear system is tri-diagonal and can be solved in O(n) using the Thomas algorithm. I am trying to solve a finite ...
4
votes
1answer
94 views

Compute specific eigenvalues in the complex plane with Feast?

In physical problems, it's quite common that we need to solve for specific eigenvalues in the complex plane, e.g. with a positive real part and negative imaginary part. In this case, we are looking ...
2
votes
0answers
52 views

Under what circumstances does Elemental's distributed SVD not work? [closed]

I am playing around with Elemental's distributed singular value decomposition and am running into two particular issues. Building the test at tests/lapack_like/SVD.cpp, and running with ...
1
vote
1answer
162 views

Appropriate Lapack/MKL routines to efficiently compute C = A* inv(B)

I know that, from numerical point of view, computing Ax = b B=inv(A), x= B*b are completely different things, and we should factor the matrix using TRF routine ...
11
votes
2answers
1k views

Compute all eigenvalues of a very big and very sparse adjacency matrix

I have two graphs with nearly n~100000 nodes each. In both graphs, each node is connected to exactly 3 other nodes so the adjacency matrix is symmetric and very sparse. The hard part is I need all ...
2
votes
1answer
59 views

Doubt regarding principled approach towards approximating the Hessian

In my optimization problem, the hessian has a structure such that it can be written as the sum of two matrices. Populating the first of the matrices is efficient. Populating the second one is ...
1
vote
0answers
61 views

What's the optimal method to solve for the top eigenvectors of a very large, real, symmetric matrix of limited rank?

Consider a real symmetric matrix of dimension N~10^5 and rank m~2000. What is the most efficient algorithm for determining the top m eigenvectors? If the answer isn't obvious, are there existing ...
9
votes
1answer
492 views

Matrix Balancing Algorithm

I have been writing a control system toolbox from scratch and purely in Python3 (shameless plug : harold ). From my past research, I have always complaints about ...
7
votes
1answer
378 views

Efficient algorithm for solving linear system with symmetric near-tridiagonal matrix?

I would like to solve the linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$, with $$\mathbf{A}=\mathbf{T}+\mathbf{C}$$ where $\mathbf{T}$ is a symmetric tridiagonal matrix and $\mathbf{C}$ is a corner-...
2
votes
2answers
611 views

Compute all eigenvectors and eigenvalues of small symmetric matrices

My problem is to compute eigenvectors and eigenvalues of a lot of small (n < 30) symetric, positive definite matrices. So far I am using LAPACK's DSYEV. The priority is speed more than accuracy. ...
6
votes
1answer
153 views

Tikhonov (Ridge) Regression and Normalization

For a typical Ridge Regression method for solving an inverse problem $$ \min_x ||A~x - b||^2 + \lambda^2||\Gamma~x||^2 $$ Which has an analytical solution of $$ \hat{x}_{est}=(A^TA+\lambda^2 \Gamma^T\...
1
vote
0answers
56 views

bound error for iterative method for solving linear system

$A$ is square and positive definite, and let $r_k = Ax_k - b$. Also let $M = \frac{1}{2}(A+A^T)$. I want to show that $$\frac{||r_{k+1}||_2}{||r_k||_2} \le \left(1-\frac{\lambda_\min(M)^2}{\lambda_\...
2
votes
1answer
95 views

linear programming feasiblity checking

Is there any sufficient and necessary condition to check the feasibility of the linear constraints $Ax=b, x\geq 0$ without solving an LP with a constant objective function? $x$ is the variable and $...
0
votes
1answer
71 views

Efficiently approximating sum of 2-norms

Suppose I have a real vector $\bf x$ of length $KN$, where $K<<N$. Let's say I break this vector $\bf x$ up into $N$ vectors each of length $K$, $${\bf x}_1,...{\bf x}_n.$$ Now I would like ...
1
vote
0answers
59 views

Efficiently computing the properties of a Chebyshev series

Suppose we have some function $f(x)$ defined as a Chebyshev series up to order $M$: $$ f(x) = \sum_{n=0}^{M} c_n T_n(x). $$ For a given coefficients vector $\mathbf{c}$, and $x \in [-1,1]$ I'm ...
1
vote
0answers
89 views

How to get the eigenvalues of Hamiltonian in an over complete basis

Let $|\psi_i\rangle$, $i=1...N+m$, be a set of overcomplete basis vector in a $N$-dim Hilbert space. The following are known: (Einstein's summation convention assumed) $$\hat{H}|\psi_i\rangle=H_{ji}|\...
1
vote
1answer
616 views

Eigen - Solving Complex Generalized Eigenvalue Problem

I've been using the Eigen C++ linear algebra library to solve various eigenvalue problems with complex matrices. I've recently had to use a generalized eigenvalue solution process, only to be ...
1
vote
0answers
116 views

Updating factorization of Laplacian (add/remove edges)

For a graph $G=(V,E)$, recall that the unweighted Laplacian is $L:=D^\top D$, where $D\in\{-1,0,1\}^{|E|\times|V|}$ is the graph "gradient" operator that subtracts adjacent vertex values onto edges. ...
2
votes
1answer
88 views

Solving system of related equations without completely recomputing LU decomposition for each equation

Let $\Sigma$ be positive definite and the $D_i$ positive diagonal. Let the $X_i$ be unknown square matrices. Consider the system of equations: $$(I+\Sigma D_i)X_i=\Sigma\hspace{5mm}\text{for}\...
8
votes
1answer
188 views

Solving two inverse problems with same solution

I've got two inverse problems, $$A_1 ~ x = b_1 \qquad A_2 ~ x = b_2$$ So far I've been solving them independently using Tikhonov Regularization and getting two estimates for $x$. However in my case $...
2
votes
1answer
366 views

Rectangular Cholesky Decomposition

Assume we have a positive semidefinite matrix $A$ with $A\in\mathbb{R^{n\times n}}$. Clearly a cholesky decomposition $A=B^TB$ exist with $B\in\mathbb{R^{n\times n}}$. For my research it would be ...
0
votes
1answer
176 views

Open Source Linear Algebra Library

I am making a code in C that requires the equivalent of Matlab's '\' command for a linear system of the form AX=B where A is an NxN matrix and X, B are Nx1 vectors- i.e a code that performs X=A\B that ...
3
votes
2answers
459 views

Python environments for AMG and Gauss Seidel as solvers instead of preconditioners

I am working on block preconditioning and seemingly it is common to write customised Krylov solvers for them. Within each solver, the individual block linear system with preconditioners are ...
2
votes
0answers
97 views

Eigenvalue decomposition of the sum: $AA^T$ + diag($u$)

Suppose $A\in\mathbb{R}^{n\times c}$,$u\in\mathbb{R}^n$,$n\gg c$. The time complexity of eigenvalue decomposing directly for matrix $AA^T+\text{diag}(u)$ is $O(n^3)$. And it is easy to avoid $O(n^3)$ ...
1
vote
0answers
30 views

Difference between Chebyshev first and second degree iterative methods

Consider linear equation $Au = f$. We want to solve it with iterative method (assuming $A$ is good). First order iterative method is: $$ u^{k+1} = u^k - \alpha_{k+1}(Au^k - f), $$ The second degree ...
7
votes
2answers
145 views

Lanczos algorithms for Hermitian system with Toeplitz kernel

Basically, I am trying to compute the SVD of a large Hermitian matrix $H$ using Lanczos iteration, while $H$ consists if a Toeplitz kernel $K$, which should be able to help speed up the matrix-vector ...
3
votes
3answers
1k views

Understanding Finite-Element Modal Analysis

I am teaching a basic course on computational physics and for the last part of the course I will introduce freshman physics undergraduates to finite-element modelling methods. I am preparing a COMSOL ...
1
vote
1answer
124 views

How to determine the truncation error with products and quotients

If I have an equation given by $$\displaystyle Y = \frac{a^2}{d^2}\frac{(1-c^2\frac{c}{a})}{(1-b^2)}$$ and I expand $a,b,c,d$ in a Taylor series, where $a$ is truncated at the $A^{th}$ order, $b$ is ...
5
votes
0answers
128 views

Rank-one Update to a Rank Revealing QR (RRQR) Factorization?

Suppose we are given an RRQR factorization for some matrix $A \in \mathbb{R}^{m \times n}$, $A\Pi = QR$ where $m > n$. Is there a cheap way to update $A' = A + uv^{\top}$ given this factorization?...
3
votes
3answers
167 views

Eigenvectors of Black-box matrix

$\DeclareMathOperator{\diag}{diag}$ Consider the generalized eigenproblem $A\mathbf{x}=\lambda B\mathbf{x}$. When solving PDEs numerically (specifically, I am interested on finding the Dirichlet ...
2
votes
2answers
516 views

Zero Eigenvalues in Lanczos Algorithm

I need to find the smallest few eigenvalues of a Hamiltonian (exact diagonalization) I use Python, and SciPy's built-in sparse eigenvalue solver. I notice, however, that for my small system (only a ...
0
votes
1answer
66 views

Show the symmetric Gauss-Seidel converges for any $x_0$

Let $A\in\mathbb{R}^{n\times n}$ is symmetric positive definite and consider solving linear system $Ax = b$. Show that the symmetric Gauss-Seidel iteration converges for any $x_0$. Solution - Since $...
0
votes
1answer
74 views

Linear stationary iteration method

Suppose you are attempting to solve $Ax = b$ using linear stationary iteration method defined by $$x_k = G x_{k-1} + f$$ that is consistent with $Ax = b$, i.e., for which $f = (I - G)A^{-1}b$. Suppose ...
2
votes
1answer
175 views

Efficient algorithm for a matrix product

Recall that a unit lower triangular matrix $L\in\mathbb{R}^{n\times n}$ is a lower triangular matrix with diagonal elements $e_i^{T}L e_i = \lambda_{ii} = 1$. An elementary unit lower triangular ...
2
votes
1answer
220 views

Determine a sufficient condition for a Hessenberg matrix to be nonsingular

Consider $A\in\mathbb{R}^{n\times n}$ whose nonzero elements are restricted to the main diagonal the strict upper triangular part, and the first subdiagonal. For $n = 8$ the locations that must be ...