Questions tagged [matrix-factorization]
Decomposition of a matrix into a product of matrices with special properties. Common matrix factorizations include LU, QR, SVD, and Cholesky.
98
questions
3
votes
0answers
103 views
Invert a huge sparse operator;
please help me with this question, I want to invert a huge sparse (non-circulant) this below in a $Ax=y$ equation:
$$(\lambda I+ \beta D+ \sigma C)x=y$$
where
I is an Identity Matrix,D is a Diagonal ...
1
vote
0answers
51 views
Using the lower triangular portion of a matrix to return a Symmetric Positive Definite Matrix?
I've recently posted this question.
To summarise, I'm dealing with supposedly Symmetric Positive Definite(SPD) matrices, but due to machine-precision they end up not being SPD.
In a comment, a user ...
3
votes
2answers
107 views
Efficient schemes for solving the extended Saddle point problem
I am interested in knowing some efficient techniques for solving the following extended Saddle point problem.
\begin{align}
\begin{bmatrix}
A & B^T & C^T \\
B & 0 & 0 \\
C & ...
4
votes
2answers
143 views
Is there a simple way to add a sparse matrix to an LU decomposition of a dense matrix?
I am solving a parabolic equation in the form:
$$
\left( {M \over\tau_j} + A \right) u^{j+1} = M f^j + {u^j \over \tau_j},
$$
where $A$ and $f$ are a dense stiffness matrix and the right hand side of ...
8
votes
3answers
524 views
Accurate Way to Calculate Matrix Powers and Matrix Exponential for Sparse Positive Semidefinite Matrices
I do need to numerically calculate the following forms for any $x\in\mathbb{R}^n$, possibly in python:
$x^T M^k x$, where $M\in\mathbb{R^{n\times n}}$ is a PSD sparse matrix, $n$ can be quite large ...
1
vote
0answers
50 views
Nonsymmetric permutations for LU factorisation of symmetric matrix
Let $A$ be a symmetric matrix. It is then well known that computing the LU factorisation of $PAP^T$ instead of $A$ for a suitably chosen permutation matrix $P$ can greatly reduce fill-in. My question ...
0
votes
1answer
68 views
Solve for large array of PD matrices
I have N matrices that are positive definite, and I have to solve for a M vectors.
As M is large in my case, doing all solves simultaneously using np.linalg.solve ...
3
votes
1answer
125 views
Quick way to find a common basis of eigenvectors between 2 matrices : valid or not?
Following the advise of @Federico Polonion a previous post, one suggested, to find a basis of common eigen vectors between 2 matrices, to simply do :
Generate 2 ...
2
votes
1answer
111 views
Efficiently compute a projection matrix from Householders reflectors
Let $A \in \mathbb{R}^{m \times n}$ where $m \geq n$.
Let $B$ and $\tau$ be the result of applying LAPACK's dgeqrfp method (R on the upper right triangle, and the ...
1
vote
2answers
130 views
Solving a specific sparse linear system without dense materialization
I need to (computationally) solve a system of equations, for the purposes of an interior point method, of the form
$$
\left[\begin{array}{cc}B & A^T \\ A & 0\end{array}\right] \left[\begin{...
2
votes
1answer
111 views
Numerical Linear Algebra: When to use Direct methods versus iterative methods to solve a linear system - for PDEs in particular
I am reading the Chapra and Canale book on numerical methods, and was working through the chapters on solving linear systems. Now the book goes through direct methods including Gaussian Elimination, ...
2
votes
0answers
65 views
Does incomplete LU preconditioning improve the asymptotic scaling of Krylov subspace methods?
It is well known that unpreconditioned Krylov subspace methods applied to the finite-difference-discretised Poisson equation with $n$ grid points per direction require $O(n \, |\log(\varepsilon)|)$ ...
1
vote
1answer
105 views
Range of a matrix from its complete orthogonal decomposition
In this StackOverflow answer, @Gokul has shown how to get a basis of the kernel of a matrix with the help of the 'Eigen' function CompleteOrthogonalDecomposition. ...
4
votes
0answers
58 views
Efficient computation of marginalized multivariate normal likelihood
In general,if we know that the marginal Gaussian distribution for some variable $\textbf{x}$ and a conditional Gaussian distribution for some $\textbf{y}|\textbf{x}$ of the forms:
$$p(\textbf{x}) = \...
0
votes
2answers
153 views
“WY” representation of QR factorization — implementations?
I have a matrix $A \in \mathbb{R}^{m \times n}$ where $m \gg n$ and I want to compute the full QR decomposition $A = QR$. Where $Q$ is an orthogonal $m \times m$ matrix.
Bishof & Van Loan (1987) ...
2
votes
1answer
128 views
Methods to improve the efficiency and the memory requirement of LU factorization for complex symmetric system matrix
I want to solve a linear set of equations (Ax=b) using LU decomposition. My "A" matrix is a complex matrix which is ...
2
votes
0answers
105 views
Computing Singular Value Decomposition of small ($4\times 4$) matrices
I need to compute the Singular Value Decomposition (SVD) of many $4 \times 4$ matrices. I'm looking for SVD algorithms specialized for small matrices. I've read that the ...
1
vote
0answers
49 views
Binary data clustering by Matrix factorization [closed]
I have read an article talking about binary clustering using Matrix factorization(see attached), but i would like to understand some optimization concepts:
Is it reasonable to use a Frobenius norm in ...
1
vote
1answer
174 views
Implementation of sparse matrix SVD for GPU
I have a sparse matrix $W$ which is almost-squared ($N+1 \times N$) and I would like to know the eigenvalues of $A = W^T W$. $A$ is Hermitian so the eigenvalues are real-positive valued.
The usual ...
1
vote
1answer
154 views
Efficient way to solve a set of linear equations $Ax=b$ when $A$ is sparse and some elements of $b$ are equal to zero
I have a set of linear equations, $Ax=b$. And about half of the elements in the right-hand side (vector $b$) are equal to zero. My system matrix $A$ is a sparse complex matrix. And $A$ is in the size ...
2
votes
1answer
79 views
Functions from Scipy, Blas, or Lapack that compute only upper triangular matrix
My goal is to transform a matrix into upper triangular form in Python. I know the function scipy.linalg.lu will do LU decomposition and get both upper and lower ...
3
votes
1answer
386 views
Why is 'scipy.sparse.linalg.spilu' less efficient than 'scipy.linalg.lu' for sparse matrix?
I have a matrix B which is sparse and try to utilize a function scipy.sparse.linalg.spilu specialized for sparse matrix to ...
5
votes
1answer
164 views
Algorithm to factorize matrix whose many rows are already of upper triangular form?
I have a matrix whose many rows are already in the upper triangular form.
$$\begin{bmatrix}
x_{11} & x_{12} & x_{13} & x_{14} & x_{5} \\
0 & x_{22} & x_{23} & x_{...
1
vote
0answers
45 views
Triangle on top of diagonal least squares
I need to solve many least squares problems with the following matrices:
$$
\pmatrix{ R \\ D_i }
$$
where $R$ is upper triangular and $D_i$ is diagonal. $R$ is the same for all the problems, while $...
0
votes
2answers
153 views
Reorder eigenvalues in Schur factorization in descending order
In this command:
[US,TS] = ordschur(U,T,select)
what should replace the select to rearrange the eigenvalues in descending ...
1
vote
1answer
98 views
How can I get Cholesky decomposition from eigenvalue decomposition?
I have
$$S = QLQ^T$$
I know $Q$, $L$, $Q^T$.
How can I get the $R$ and $R^T$ for the Cholesky decomposition $S=R^TR$?
2
votes
1answer
158 views
Bareiss algorithm vs. LU-decomposition
I at the moment try to fully understand the Bareiss algorithm for calculating determinants.
One question that came to my mind is the following:
Why is LU-decomposition much more often used than the ...
5
votes
0answers
65 views
Why the two Gram-Schmidt algorithms produce different results for qr factorization?
For the qr factorization using classic Gram-Schmidt algorithm, I found the 2 different implementations below. The first one uses the for loop to compute the upper ...
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
2answers
224 views
Block-matrix: optimal fill-in reduction for LU factorization
Consider a square $N \times N$ block-matrix $\mathbf{A}$, where each $n \times n$ block $\mathbf{A}_{ii}$ is either a dense block or a zero-block. So, $N$ denotes the number of blocks, $n$ denotes the ...
2
votes
1answer
205 views
Reference for QR algorithm for complex matrix
I am trying to find out if the known QR algorithm to find the eigenvalues of a real matrix, which can be found in the book Fundamentals of Matrix Computations, can also be used for complex matrices ...
0
votes
0answers
88 views
Why does the matlab command **chol(A)** slower than **chol(A,'lower')** for a large sparse SPD matrix?
For a SPD matrix A, there exists Cholesky factorization $A=LL^T$ or $A=R^TR$, where L, R are a lower and upper triangular matrix, respectively.
Also in matlab, there has a command R = chol(A) which ...
7
votes
3answers
416 views
Does a symmetric positive definite matrix also have $\mathbf{A} = \mathbf{L}^T\mathbf{L}$ (where $\mathbf{L}$ is a lower triangular matrix)?
As we know, for a symmetric positive definite (SPD) matrix $\mathbf{A}$, there is a theorem about the Cholesky factorization $\mathbf{A}= \mathbf{L}\mathbf{L}^T$, where $\mathbf{L}$ is a lower ...
4
votes
1answer
148 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}$
...
4
votes
1answer
148 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
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
125 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
1answer
692 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 ...
1
vote
1answer
177 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
208 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
1answer
133 views
Classical vs. modified Gram-Schmidt
It is often said that modified Gram-Schmidt is more robust with respect to rounding errors than classical Gram-Schmidt, but it is very hard to find a good explanation / example of why this is so. Can ...
3
votes
1answer
426 views
Givens rotation vs 2x2 Householder reflection
The usual story of Givens rotations vs Householder reflections is that Householder reflections are better if you want to map a long vector to $e_1$, while Givens is better if you want to map a 2-...
1
vote
1answer
104 views
Software for parallel incomplete LU factorisation
I am looking for a software package to compute incomplete LU factorisations in parallel. Further considerations are:
The package must allow for arbitrary level-of-fill or threshold-based truncation. ...
2
votes
0answers
223 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|\...
4
votes
2answers
260 views
Inverting a matrix from LU decomposition
The LAPACK routines xGETRI compute the inverse of a matrix $A = PLU$ in its LU decomposed form by first computing $U^{-1}$, and then solving the system:
$$
(A^{-1} P) L = U^{-1}
$$
My question is: ...
2
votes
1answer
265 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 ...
1
vote
1answer
221 views
Pivoted Cholesky vs Modified Cholesky
I am solving nonlinear least squares problems with the normal equations approach, so on each iteration, I need to solve:
$$
J^T J \delta = -J^T f
$$
for the step $\delta$, where $J$ is a large (...
2
votes
0answers
42 views
ILUTP in sparse.linalg.spilu?
In Matlab, an ILU with threshold and pivoting (ILUTP) can be passed by default as:
setup.type = 'ilutp';
[L, U] = ilu(A, setup);
Looking for an equivalent in ...
3
votes
0answers
72 views
Numerical analysis, pivoting and incomplete LU decomposition
When doing LU decomposition, the algorithm will break down if any of the diagonal element $x_{ii}$ is zero. Therefore, we can use pivoting on the matrix such that $x_{ii}$ is no longer zero. That is ...
4
votes
1answer
620 views
Numerically find the nearest positive semi definite matrix to a symmetric matrix
I have a symmetric matrix $M$ which I want to numerically project onto the positive semi definite cone.
To do so, I decompose it into $M = QDQ^T$ and transform all negative eigenvalues to zero. (...
