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.
976
questions
0
votes
1answer
63 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
119 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
89 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
124 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
66 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
209 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
109 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
95 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
61 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
188 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
91 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
161 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
103 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,...
2
votes
1answer
63 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 ...
4
votes
1answer
105 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
94 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$ ...
3
votes
2answers
536 views
Fast and accurate eigenvalue computation for 3x3 posdef matrices
I am looking for a very fast and efficient algorithm for the computation of the eigenvalues of a $3\times 3$ symmetric positive definite matrix.
The algorithm will be part of a massive computational ...
4
votes
2answers
145 views
Singular values of $X$ in $AX+XA=C$?
Suppose I have semi-positive definite matrices $A$ and $C$, is there an efficient approach to get top singular values of X entering the following expression?
$$
AX+XA=C
$$
My matrices are 4k-by-4k ...
4
votes
1answer
341 views
Diagonalize a unitary matrix with orthogonal matrices using numpy
An important component of the Cartan KAK decomposition for 2 qubit operations is to diagonalize a 4x4 unitary matrix using orthogonal (not unitary, purely real orthogonal) matrices. That is to say, ...
1
vote
2answers
327 views
Fastest algorithm for pseudoinverse of skinny matrices
For a performance-sensitive problem, I need to compute the pseudoinverse of a skinny matrix (#rows = 1000–10000, #cols= 10–20).
I already employ the traditional SVD econ method. For some problem ...
2
votes
1answer
347 views
Implementation of Jacobi iteration
I have implemented the Jacobi iteration in C++ using a dense vector and a sparse matrix in CSR format. The code is as follows:
...
2
votes
1answer
473 views
Cholesky for ill-conditioned/singular covariance matrices
Can someone suggest a way to get Cholesky factorization of a singular covariance matrix? I need it to match Cholesky on full-rank matrices, ie coordinate order should be preserved. My attempt below ...
2
votes
0answers
41 views
Inverses of many standard subspaces of one large matrix
i have a large rectangular invertible matrix M (about 5000x5000) and i have a loop in which i do the following for each iteration i (there are about 6000 iterations):
i am given a subspace S_i (which ...
7
votes
0answers
118 views
Implementation of Lanczos method that returns tridiagonal matrix
The Lanczos method can be used to obtain extremal eigenpairs of sparse symmetric or hermitian matrices. I know there are several implementations of the Lanczos method (as well as Arnoldi, Davidson, ...
1
vote
1answer
117 views
Incomplete LU decomposition of sparse matrix
I have a sparse matrix stored in CSR format. For this matrix, I would like to get the incomplete LU decomposition. I tried to find algorithms which can utilize the CSR format but I could not find ...
1
vote
1answer
145 views
How to use QZ decomposition for single matrix in Matlab?
Can I use QZ decomposition on a single square matrix in Matlab?
Like,
[Aa,Q,Z]=qz(A);
1
vote
2answers
95 views
Fast iterative approximate order-oblivious Orthogonalization algorithm?
I have set of N m-dimensional vectors $\{\phi_i\}$ which gradually loose mutual orthogonality in an algorithm. => I have to re-orthogonalize them every few iterations. But if I do e.g. Gram–Schmidt ...
3
votes
1answer
75 views
Why are the round-off errors when solving the linear system $Ax = b$ of order $\varepsilon_\text{mach} x_j$?
I was reading a paper on arXiv where, in Section 2.4, the authors are discussing the error that arises in the solution of a linear system
$$Ax = b,$$
or, to match up better with the paper,
$$\Phi \...
4
votes
1answer
403 views
Is there a library that allows einstein summation on dense, sparse, and LinearOperator type tensors
Numpy's einsum only works with dense tensors.
Is there an alternative that also works with sparse tensors and linear operators?
For example, I might have a ...
3
votes
1answer
161 views
Matrix condition number and reordering
Does the condition number change when a matrix is reordered by e.g. Cuthill Mckee or some other method?
5
votes
0answers
253 views
Is there any catch on using `zgemm3m` vs regular `zgemm`?
I've just (to my embarrassment) encountered a BLAS-like extension of a matrix-matrix product subroutine gemm in Intel MKL: gemm3m...
1
vote
1answer
92 views
How to implement this even-odd matrix decomposition efficiently?
Note: This question has also been asked on stackoverflow - see https://stackoverflow.com/questions/57197910/how-to-implement-this-even-odd-matrix-decomposition-efficiently?noredirect=1#...
1
vote
0answers
48 views
Does Boost provide a template implementation of the wedge product?
Does the boost C++ library implement the computation of the wedge product?
The wedge product is mentioned here, but it is not very clear (to me at least) whether there is a template implementation of ...
1
vote
0answers
81 views
What's the more efficient way to solve this matrix equation?
This is intended to be a more generic question not about a specific system. Given a hermitian matrix $H(x_1,\dots,x_n)$ depending non-linearly on some real parameters $x_1,\dots,x_n$. We want these to ...
2
votes
1answer
56 views
How to pass matrices to parallel workers quickly in matlab?
I am trying to solve many different linear systems in parallel in matlab. The problem is, each linear system has entirely different parts and are fairly large, so passing the information to each of ...
1
vote
1answer
119 views
Sparse matrix inversion
I have the impedance matrix $Y$, formulated from an electrical network by augmented nodal analysis. The matrix $Y$ is shown as an image to illustrate its feature visually, where all the white blocks ...
2
votes
1answer
60 views
numpy.outer without flatten
$x$ is an $N \times M$ matrix.
$y$ is a $1 \times L$ vector.
I want to return "outer product" between $x$ and $y$, let's call it $z$.
z[n,m,l] = x[n,m] * y[l]
...
0
votes
1answer
83 views
what does “D = diag(W.1)” means?
, what does “D = diag(W.1)” means?on page #2, just below equation (6)
PFA screenshot and here is the link of the paper -
original paper
5
votes
1answer
350 views
Fast algorithm for computing cofactor matrix
I wonder if there is a fast algorithm, say ($\mathcal O(n^3)$) for computing the cofactor matrix (or conjugate matrix) of an $N\times N$ square matrix. And yes, one could first compute its determinant ...
2
votes
1answer
113 views
Null space for smoothed aggregation algebraic multigrid
I do not really get the point of null space usage for creating the prolongation operator for smoothed aggregation algebraic multigrid. I know what the null space is per definition and I know that the ...
2
votes
0answers
164 views
Regularized least squares with QR factorization
Consider the regularized least squares problem
$$
\min_x || b - A x ||^2 + \lambda^2 ||x||^2
$$
which is equivalent to
$$
\min_x \left|\left| \pmatrix{b \\ 0} - \pmatrix{A \\ \lambda I} x \right|\...
5
votes
1answer
147 views
Complex differentiation of linear solvers
I have a linear system $$Ax=b$$ which I'm solving approximately, and I need to take the frechet derivative of x with respect to z. Were I solving the problem exactly (either analytically or to machine ...
3
votes
1answer
179 views
Optimality of block-Jacobi preconditioner
For a dense $N \times N$ matrix $A$, is the block-Jacobi preconditioner comprising the inverse of the diagonal blocks of $A$ the optimal block-diagonal preconditioner? Could there exist another matrix ...
1
vote
0answers
56 views
Kinetic preconditioning
Publication arXiv:0804.2583 describes a method for doing self-consistent iteration without having to diagonalize the Hamiltonian operator at every step.
IX. PRECONDITIONING
As already mentioned, ...
0
votes
0answers
34 views
A preconditioner for self-consistent iteration
I tried to derive a preconditioner for self-consistent iteration similar
to section IX in arXiv:0804.2583.
For simplicity, consider here only
one orbital (one or two electrons) systems.
Suppose that ...
3
votes
1answer
128 views
Whitening transformation does NOT return a unit covariance matrix
For this question, I am using the following Wiki definition of Matrix whitening:
Suppose $X$ is a random (column) vector with non-singular covariance matrix $\Sigma$ and mean 0. Then the ...
2
votes
1answer
149 views
In iterative methods, are matrix decompositions considered useful for implementation?
When we study an iterative method from textbooks, for example, see the Gauss-Seidel Method, the given matrix is decomposed with suitable splittings. In the example, $A = L+U$. So we can proceed with ...
5
votes
2answers
398 views
Python: vectorizing a structured linear system solve
Overview
I am looking for a way to solve a structured linear system in Python without using a for loop (preferably using vectorization, if possible).
Background
Consider the following linear system:...
2
votes
1answer
185 views
Re-using LU factorization within iterative (?) setup for a sum of two matrices
So, I would love to make at least some use of my preexisting data, no matter how small, and just out of ideas. Maybe I am just a prisoner of a Kahneman-like theatre-ticket paradox, and don't know ...
4
votes
1answer
251 views
Accurate way of getting the square root inverse of a positive definite symmetric matrix
What is the most accurate algorithm to get the square root inverse of a positive definite symmetric matrix? I am not looking as much for efficiency, though using quadruple precision computation is out ...
