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
9 votes
1 answer
197 views

What makes a good preconditioner when only a few approximate iterations are needed?

For deterministic solver of $Xw=y$, one recommendation is to pick $P$ such that $P^{-1}X$ has a low condition number. However, this condition only really matters when you want to reduce initial error ...
3 votes
2 answers
735 views

Why does this preconditioner effectively reduce the condition number of a random SPD matrix?

Consider some randomly generated matrix $B\in\mathbb{R}^{100\times100}$ and let $A:=BB^{\top}$ On MATLAB I computed the condition number of $A$, I obtained a value of $2.8377\mathrm{e}+04$ However if ...
1 vote
0 answers
63 views

Does cblas_dgemm mutate my input matrices?

I have written a matrix class Matrix<T> for which I have implemented a wrapper function for cblas_dgemm. ...
2 votes
0 answers
112 views

Global minimum of a function involving matrix exponentials

I have a matrix system \begin{equation} \vec{x}(t) = e^{(iAt)} \vec{x}(0) \end{equation} that comes from the solution of the differential system \begin{equation} \frac{d}{dt}\vec{x} = iA\vec{x}, \end{...
40 votes
23 answers
10k views

Good examples of "two is easy, three is hard" in computational sciences

I recently encountered a formulation of the meta-phenomenon: "two is easy, three is hard" (phrased this way by Federico Poloni), which can be described, as follows: When a certain problem is ...
8 votes
2 answers
231 views

Is there one general approach to build a projection methods for different problems?

My question is probably going to be too general to answer it with a couple words. Could you please suggest a good reading in that case. Projection methods are used to reduce size of the solution space ...
11 votes
3 answers
507 views

Parallel algorithm for eigensystem of a tridiagonal matrix

I'm doing a Lanczos diagonalization of a large sparse matrix (~2 million elements). Almost all of the steps in the Lanzcos algorithm are done in parallel on the GPU, except for diagonalizing the ...
3 votes
0 answers
54 views

Sparse generalized symmetric eigensystem solver

Can anyone recommend a good software for solving generalized symmetric eigenvalue problems of the form, $$ A x = \lambda B x $$ where $A,B$ are symmetric and sparse, and $B$ is positive definite? I ...
2 votes
0 answers
48 views

Finding spectrum of a Kronecker factored + block-partitioned matrix

I have dense $d\times d$ matrices $A$, $B$, $C$ with $d\approx 1000$ and want to find the top $10^5$ eigenvalues of the following positive definite matrix: $$ \Sigma= \left(\begin{matrix} A\otimes A &...
4 votes
2 answers
116 views

Numerical estimation of eigenfunctions of Laplacian

Consider the Laplace equation, $$ \nabla^2 f(r,\theta,\phi) = 0 $$ in spherical coordinates. We know that the solution to this equation can be derived analytically, and is given by, $$ f(r,\theta,\phi)...
6 votes
1 answer
201 views

Which preconditioners make Richardson iteration convergent?

Suppose we solve an $m\times n$ full-rank system of equations $Ax=b$ by iterating the following for a small enough $\mu>0$ $$x=x+\mu B(b-Ax)$$ Is there a nice description of kinds of $B$ which make ...
2 votes
1 answer
140 views

Matrix regularisation for ill-conditioned problems

I have read that matrix regularisation can improve the stability of LU or Cholesky decomposition of ill conditioned problems. The idea is to add a value to the diagonals of a matrix: $B=A+cI$ In the ...
3 votes
1 answer
138 views

Column-normalized inverse?

Suppose we define $A^{*}$ of positive definite $A=X'X$ using following two steps: let $B=A^{-1}$ scale columns of $B$ to obtain a matrix with $1$'s on the diagonal For the case of singular $A$, we ...
7 votes
0 answers
61 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^\...
5 votes
1 answer
101 views

Largest singular value without using the adjoint

The square of the largest singular value of a linear map $A$ can be computed by using the power iteration for $A^TA$ and one advantage of this is that the iteration is matrix free, i.e. you only need ...
2 votes
1 answer
114 views

Numerical representation of linear spaces

A linear space in mathematics is a set whose elements (vectors) have operations of addition and multiplication by a scalar defined in such a way that certain properties are satisfied (commutativity, ...
3 votes
1 answer
1k views

Set of linear ordinary differential equations with a mass matrix

What methods are known for efficiently solving a large set of linear homogeneous ordinary differential equations of the following form? \begin{equation} \mathbf{B} \frac{d\mathbf{y}}{dt} = \mathbf{A} \...
9 votes
2 answers
2k views

Is there an algorithm or graph theory that allows me to not need to store an intermediate matrix when calculating AT*Y1*A + BT*Y2*B?

I have a system of conductors for which there are two dense matrices of the (complex) mutual admittances, $Y_A$ and $Y_B$, which are symmetric. Then, an equivalent nodal admittance matrix $Y_N$ is ...
7 votes
3 answers
2k views

Choosing subset of vectors to approximate a subspace

Suppose I have a high-dimensional vector space $X$, a subspace $V \subset X$, and a collection of $n$ vectors $\{x_i\}_{i=1}^n \subset X$. My question is: How can I choose a small collection $k < ...
3 votes
1 answer
168 views

Choice of iterative solver for a sparse asymmetric matrix with symmetric structure

I have a sparse $nxn$ matrix A with pretty interesting structure. It has a block structure with symmetric structure but asymmetric blocks. Expressed mathematically the block $A_{jk} = A_{kj}$ but $A_{...
57 votes
4 answers
9k views

What guidelines should I follow when choosing a sparse linear system solver?

Sparse linear systems turn up with increasing frequency in applications. One has a lot of routines to choose from for solving these systems. At the highest level, there is a watershed between direct (...
3 votes
1 answer
85 views

Using submatrices of matrix decomposition for solving a large number least-squares problems

I want to decrease the computational time for solving a large number (>1000) of least-squares problems. Given a matrix, the system matrix for each least-squares problem is a submatrix of the given ...
1 vote
1 answer
94 views

constructing a symmetric matrix for finite difference

I come across the following operator in a paper $\mathcal{I}\psi = \psi_{xxxx} + (r~\psi_x)_x$, where $\psi=\psi(x)$ and $r=r(x)$. Periodic boundary condition is employed. It claims that the operator $...
10 votes
2 answers
567 views

Extracting diagonal of an approximately diagonal matrix when we don't know its entries

What is a good way to extract the diagonal from a symmetric matrix that is already almost diagonal, but where you don't have the matrix elements (only the ability to apply it to vectors)? Further ...
4 votes
1 answer
179 views

Efficient Triangularisation of $\mathbf{S} = \operatorname{triag}\left(\mathbf{A}\right)$ ; $\mathbf{S}\mathbf{S}^T = \mathbf{A} \mathbf{A}^T$

The most computationally intensive part of my application is the triangularisation of a matrix, $\mathbf{S} = \operatorname{triag}\left(\mathbf{A}\right)$, such that $\mathbf{S}\mathbf{S}^T = \mathbf{...
0 votes
0 answers
37 views

Solving Symmetric/Hollow Matrix issue

I have a particular issue and need something creative or solution from calculus. I have Symmetric/Hollow Matrix, a numbers are % of mismatch between them. Ideally, all of them should be 0, but I have ...
5 votes
2 answers
332 views

Do most statistical packages and libraries in high-level programming languages rely on LAPACK for their matrix inversion operations?

Possible an open-ended question, but I am wondering if most statistical packages and libraries, for instance, Stata, R, Python's NumPy and MATLAB rely on LAPACK algorithms to perform matrix operations,...
0 votes
0 answers
71 views

How to estimate stability and stiffness of a system of coupled ODEs?

I'm running into issues with Python/Julia ODE solvers requiring prohibitively small timesteps to evolve a system of 4 coupled ODEs (the order of magnitude of the state variables and time unit span ~40-...
1 vote
2 answers
68 views

Calculate average distance between pairs of points without computing full distance matrix

Suppose I have a set of $N$ points of shape $N \times D$, where $D$ is the dimensionality. I want to compute the average Euclidean distance between all points, as well as additional moments such as ...
2 votes
1 answer
92 views

Finding the correct order of eigenvectors of a parameter-dependent Hermitian matrix

so, I have a symmetric, analytic matrix $\mathbf{H}(x)$ ($x$ is real). Because $\mathbf{H}(x)$ is analytic and $x$ is real, it is possible to find analytic functions for the eigenvectors and the ...
-1 votes
1 answer
39 views

Finding matrix of a linear transformation using R programing

The full question is: Let {u1, u2,···, un} and {x1, x2,···, xm} be bases for Rn and Rm respectively. Let T:Rn→Rm be the linear transformation whose associated matrix with respect to the above bases is ...
2 votes
0 answers
88 views

Closed-form Jacobian of SE(3) element with 6-degrees of freedom

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=AB^{-1}$. $W$ can also be ...
2 votes
1 answer
106 views

Numerical diagonalization of Hamiltonian

Framework I am trying to diagonalize the Bogoliubov-de Gennes Hamiltonian. The problem is that the Hamiltonian contains a Laplacian. This could be solved by using a discretized Laplacian. How I tried ...
3 votes
1 answer
63 views

Requesting for Finite Difference Methods reference in Portuguese or English

Crossposted on Mathematics SE I have been assigned a group project for an introductory Linear Algebra subject on Finite Difference Methods and sparse matrices. Our professor advised we use Gilbert ...
0 votes
0 answers
38 views

Matrix 3D rotation/translation/scaling relative to each other

I am trying to implement an algorithm, that removes a part of an assembly (3D) and translates/rotates/scales those other matrices which are child/parent to a new matrix. For example: (MT = ...
4 votes
1 answer
218 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 ...
2 votes
0 answers
43 views

Solving linear system and obtaining operator norm

I need to solve a linear system of the form $(\mathrm{Id} + \mathbf{J})\mathbf{x} = \mathbf{b}$ for $\mathbf{x}$ and I also need to compute the operator norm of $\mathbf{J}$ (i.e. the largest singular ...
1 vote
1 answer
75 views

Fast algorithm to compute chi-square

I would like to evaluate the chi-square of the form $\chi^2=v^{T}C^{-1}v$ where $v$ is a column vector and $C$ is a covariance matrix. Both $v$ and $C$ are known and $C$ is a $740\times740$ matrix. ...
0 votes
0 answers
175 views

Solving huge dense square symmetric linear system

I have a linear system of the type $A x = y$ where A is a dense, square, symmetric, positive definite matrix, $x$ a vector of unknown parameters, and $y$ is a vector of observed quantity. I know that ...
6 votes
1 answer
114 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 ...
5 votes
1 answer
190 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 ...
9 votes
2 answers
246 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 \...
16 votes
3 answers
1k views

Why would a computational scientist need to implement their own version of std::complex?

Many of the better-known C++ libraries in computational science such as Eigen, Trilinos, and deal.II use the standard C++ template header library object, ...
7 votes
0 answers
102 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 ...
0 votes
0 answers
85 views

Weighted Jacobi Not Working on 1D Poisson (Issue with Optimal $\omega$)

I've been trying to learn some numerical linear algebra, and I decided to try to implement the weighted Jacobi method to the 1D Poisson problem $$-u''(x)=f(x),\qquad u(0)=a,\ u(1)=b,$$ where we ...
2 votes
2 answers
70 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. ...
5 votes
0 answers
121 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 ...
0 votes
1 answer
50 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 ...
2 votes
2 answers
2k views

What sparse linear programming solver it is better to use?

I have the following LP problem: $$ \min \limits_{\varepsilon, x_{1}, \ldots, x_{n}}f(\varepsilon, x_{1}, \ldots, x_{n}) = \varepsilon \;\;\;\;\; \mathrm{s.t.} \;\;C x \geq 0, \;\; x_{i}^{0} - \...
1 vote
1 answer
431 views

Direct or iterative solver for ill-conditioned problems

I have to solve an ill-conditioned sparse matrix. Once I read that iterative solvers are the better tool for such problems. Is that true? If yes, why?

1
2 3 4 5
22