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,010
questions
2
votes
3answers
492 views
Generating Random Orthogonal Matrices in C++
I'm looking for an open-source library for the generation of random n-dimensional orthogonal matrices in C++.
In python, it looks like such a function is available in the NumPy package. But I was not ...
2
votes
1answer
205 views
Sparse matrix-matrix multiplication using AVX2
I have two sparse general matrices stored in CSR format I need to multiply. Is there any chance to gain performance using AVX2? In general the matrices are big (hundreds of millions of non-zeros and ...
0
votes
1answer
344 views
How to implement the Hessenberg QR Algorithm?
For context, I'm creating a linear algebra library from scratch for learning purposes in C. Right now I'm working on calculating eigenvalues but my implementation of the QR Algorithm is diverging. ...
9
votes
1answer
2k views
Increasing computational performance by using 16 bit numbers
I recently found the following article where it was stated that using 16 bit numbers can be used to increase the computational performance of AI applications. According to the article numbers above 16 ...
0
votes
2answers
227 views
Update for QR factorization least squares
I found after some research that the most numerically stable way to solve the least squares problem is through QR factorization. For $n$ number of observations and $p$ number of parameters it takes ...
3
votes
2answers
1k views
Cholesky decomposition vs LDL decomposition
In different books and on Wikipedia, you can see mentions of Cholesky decomposition and only sometimes of LDL decomposition.
As far as I understand, LDL decomposition can be applied to a broader ...
3
votes
1answer
116 views
Big Theta Complexity of Gaussian Elimination using Complete Pivoting
I already know the Big O for partial pivoting is $O(n^3)$ and remain the same for complete pivoting. I also know the big theta complexity for partial pivoting is $2/3 n^3$
I would like to know the ...
2
votes
0answers
92 views
Solving a huge least squares system of equations when I can only evaluate Ax
I have a situation where I can generate a system of $M$ linear equations for $N$ variables ($N \ne M$). Implicitly this is of the form $Ax=b$ with $A \in \mathbb{R}^{M \times N}$, although I never ...
3
votes
1answer
53 views
Nonlinear least squares resolution matrix
For a linear least squares problem, it is possible to define a resolution matrix, relating the estimated model parameters to the true model parameters. If we are solving a regularized problem,
$$
\...
5
votes
0answers
137 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 ...
1
vote
2answers
151 views
Why Householder transformation can not be chosen to be an identity matrix?
For Householder transformation, we know that
$H = I-uu^T$, where $\|u\|_2=\sqrt{2}$. When it acts on any vector $x$, $Hx$ and $x$ is symmetric with respect to $span(u)^T$. But I have read a ...
3
votes
1answer
2k views
Efficient ways to numerically evaluate matrix exponentials
What are some computationally efficient ways to solve matrix exponentials, i.e. functions of the form :
$f(X)=e^{X}$, where $X$ is a square matrix?
So far I have been able to diagonalise some ...
1
vote
1answer
43 views
Determinant of a matrix after removing or adding lines and columns
In quantum mechanics, the wavefunction of N electrons is given by a determinant. I am working on a Monte Carlo algorithm. At each Monte Carlo step, I need to add or remove an electron, which means ...
3
votes
2answers
704 views
Inverting really big symmetric block diagonal matrix
I have a really big symmetric 7.000.000 X 7.000.000 matrix that i would like to invert. The matrix is extremely sparse and it can be rearranged as to become a block diagonal matrix. The biggest blocks ...
35
votes
2answers
12k views
why is A*v+B*v faster than (A+B)*v?
$A$ and $B$ are $n \times n$ matrices and $v$ is a vector with $n$ elements. $Av$ has $\approx 2n^2$ flops and $A+B$ has $n^2$ flops. Following this logic, $(A+B)v$ should be faster than $Av+Bv$.
Yet,...
1
vote
0answers
65 views
Can Representation Theory be studied computationally / numerically?
Can a subfield such as the representation theory of Lie algebras be studied computationally / numerically -- is there an interplay between the abstract and the concrete? I would be grateful for an ...
2
votes
1answer
136 views
Solution method of nonlinear heat transfer analysis
The governing equation of transient heat transfer analysis is described as follows:
$$C \frac{dT}{dt}+K T = Q$$
When using backward difference scheme for the discretization of the time we get the ...
5
votes
0answers
92 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 ...
2
votes
0answers
55 views
Avoid matrix multiplication in algebraic multigrid method
Currently when I try to solve a linear algebra system of the form of $A x =b$ I use the algebraic multigrid method. The algebraic multigrid method uses a Galerkin product to form the coarse grid ...
5
votes
0answers
91 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 & ...
1
vote
0answers
60 views
Numerical methods. MDF (ILU) implementation
I am trying to implement Minimum Discarded Fill (MDF) Ordering algorithm for incomplete matrix factorization. The algorithm description is here on page 60 Preconditioning Techniques for a Newton–...
5
votes
1answer
375 views
Fastest way to calculate the $2$-norm (or an upper bound for the $2$-norm) of the inverse of a matrix $A\in \mathbb{C}^{N\times N}$
I have a matrix $A\in \mathbb{C}^{N\times N}$ and I need to calculate $||A^{-1}||_{2}$ efficiently. Can it be done without having to evaluate the inverse explicitly?
In general, I am looking for ...
2
votes
0answers
48 views
Fitting a plane with the Prewitt gradient operator
Prewitt gradient operator
Show that the Prewitt gradient operator can be obtained by fitting the least-squares plane through the 3 × 3 neighborhood of the intensity function.
Hint: Fit a plane to ...
10
votes
3answers
616 views
Linear algebraic research direction that's not to do with differential equations and physics?
So I've found some interesting linear algebraic research areas that's both pure-ish, with a numerical bent to it, too -- e.g. inverse eigenvalue problems have both interesting theoretical and ...
5
votes
0answers
50 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
$$...
4
votes
3answers
176 views
Solving for a vector in a linear system that is both left and right multiplied
I have a linear system where I am given 2 matrices, $A$ and $B$, and 2 vectors, $v$ and $c$, and I need to solve for the vector $x$. $A$ is $n\times n$, $B$ is $n \times n \times n$, and the vectors $...
3
votes
0answers
180 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 ...
4
votes
1answer
83 views
Is operation count a reliable predictor of performance when comparing two formulations?
I have two formulations to solve a problem (both give dense, complex and symmetric system). They are solved multiple times in a loop. I am trying to predict which is better to use.
The first one ...
3
votes
0answers
44 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 ...
1
vote
1answer
98 views
Spline regularization
I am fitting some B-splines to data, but the data has a "gap" region where the spline is less constrained by the data. I want to devise a regularization scheme to help prevent the spline from ...
0
votes
1answer
68 views
Method to calculate solution of a linear equation system?
I am searching a solution method for the following equation system of equation systems:
Let $A, B \in \mathbb{R}^{n \times n}$ be s.p.d. Matrices and $O$ be the zero matrix of the same size. Further ...
1
vote
1answer
344 views
Solve linear system with Newton-Raphson method
Is it possible to solve a linear matrix system $A x = b$ using the Newton-Raphson method? If yes, how can this be done? More special, how is the derivative build?
4
votes
1answer
153 views
Low rank update of QR of inverse
I am in a situation where as part of a sort of inverse power method scheme, I want to very often perform the following step:
Apply a symmetric rank one update $uu^\top$ to my inverse matrix $A^{-1}$
...
0
votes
1answer
43 views
How to find the nearest point inside a list in a given direction
Being $\bar{\mathbf{x}} \in \mathbb{R}^3$ a point and $S =\{\mathbf{x}\}_{i=1}^N \in \mathbb{R}^3$ a sample of N points. I am looking for a simple algorithm to determine the nearest point in $S$ in ...
0
votes
1answer
65 views
Analytic formula for $\arg\max_{\|z\|_\infty \le 1}z^T A z$, where $A=uu^T+vv^T$
Let $u$ and $v$ be column vectors of size $n \gg 1$ (not both zero), and consider the matrix $A:=uu^T+vv^T$
Question
What is an analytic formula for $\arg\max_{\|z\|_\infty \le 1}z^TAz=\arg\max_{\|z\...
3
votes
1answer
140 views
Analytic formula for leading eigenvector of $uu^T + vv^T$?
Let $u$ and $v$ be nonzero column vectors of size $n$ and consider the $n \times n$ positive-definite matrix $A:=uu^T + vv^T$. In this post https://math.stackexchange.com/a/112201/168758, the ...
3
votes
1answer
109 views
Efficient computation of leading eigenvector of a matrix product of the form $ADA^T$, where $D$ is diagonal
Let $A=[A_1|\ldots|A_m] \in \mathbb R^{n \times m}$ with $n \gg m \gg 1$ and $D=\text{diag}(d_1,\ldots,d_m)$ where $d_1,\ldots,d_m > 0$, and consider the $n\times n$ positive-definite matrix $X=\...
4
votes
1answer
157 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,...
1
vote
0answers
96 views
implementation for coppersmith matrix multiplication
Is there any online implementation for the coppersmith matrix multiplication I have searched alot but can not find any? and if there is not any why is that Isn't this algotithm much faster than ...
1
vote
0answers
313 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-...
1
vote
1answer
126 views
Lapack symmetric update $B^{-1}AB^{-T}$
Does Lapack have a routine that, given symmetric $A=A^T$ and $B$, computes the symmetric matrix $B^{-1}AB^{-T}$ (while preserving symmetry exactly)?
It would be enough to have this routine for ...
1
vote
1answer
118 views
Vectorization of Jacobi iteration
Assume I have a linear system of $A x = b$ which I want to solve with Jacobi iteration. Matrix $A$ is given in CSR format. The vectors are dense.
The code for Jacobi iteration is quite clear and can ...
4
votes
1answer
88 views
Value of $\gamma$ in the H-infinity norm
Suppose I have the system:
$$\dot{x} = Ax+Bu\\
y=Cx+Du$$
and the following Hamiltonian matrix:
$$H=\begin{pmatrix}
A & \frac{1}{2}B^TB\\
-CC^T&-A
\end{pmatrix}$$
I want to find the ...
5
votes
0answers
193 views
Solving $AXB + X\odot C = D$ matrix equation
Can anyone see a way to solve this equation efficiently?
$$AXB + X\odot C = D$$
I tried a straightforward solution that involved vectorizing $X$ but that turned out too expensive for my application -...
1
vote
1answer
103 views
Using LAPACK to compute $B^{-1}AB^{-T}$ for thin $B$
How can I use BLAS/LAPACK to compute
$$
B^{-1}AB^{-T}
$$
where $A\in\mathbb{R}^{n,n}$, $B\in\mathbb{R}^{m,n}$ is full rank matrix with $m>n$, and $B^{-1}y:=\arg \min_{x} \|Bx-y\|_{2}$.
In theory, ...
6
votes
1answer
177 views
Computing square root of diag(u)-uu'?
I need an efficient way to take square root of a matrix which is a sum of diagonal matrix and rank-1 matrix.
More specifically it's the following matrix
$$A=D-uu'=\text{diag}(u)-uu'$$
Where entries ...
2
votes
1answer
112 views
Efficiently finding binary vectors satisfying multiple conditions
I am trying to solve the following problem:
Given a binary matrix $\mathbf{A} \in \{0,1\}^{m \times n}$ and a vector $\mathbf{b} \in \mathbb N^n$, does there exist a binary vector $\mathbf{c} \in \{0,...
1
vote
1answer
68 views
How to avoid unnecessary checks when inverting this LU decomposition
Background for the question
I am currently working on a Matlab code in which the systems of linear equations $Ax_1 = b_1$, $Ax_2 = b_2$, ... have to be solved. As the matrix $A$ is constant during ...
3
votes
1answer
130 views
Weighted QR Implementation
Say I want a QR decomposition of matrix $A$, where orthogonality of $Q$ is with respect to a generic non-degenerate positive-definite bilinear form $\phi$ (in my case, $\phi$ is "defined" by a finite-...
2
votes
0answers
105 views
Small residual but wrong results
When I use BiCGStab to solve a linear matrix system, I use the relative residual to exit the iteration and output the results. For calculating the relative residual I divide the norm of vector $r$ ...