Questions on the algorithmic/computational aspects of linear algebra, including the solution of linear systems, least squares problems, eigenproblems, and other such matters.

learn more… | top users | synonyms

0
votes
1answer
51 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 ...
1
vote
1answer
46 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 ...
3
votes
0answers
65 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
23 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 ...
3
votes
1answer
35 views

Lanczos algorithms for Hermitian system with toeplitz kernel

Basically I am trying to compute the SVD of a large hermian 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
1answer
116 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 ...
0
votes
1answer
36 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 ...
3
votes
0answers
28 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
2answers
72 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 ...
3
votes
2answers
131 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
40 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
66 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 ...
0
votes
0answers
13 views

Column form elementary matrix in terms of element form elementary matrices

Recall that any unit lower triangular matrix $L\in\mathbb{R}^{n\times n}$ can be written in factored form as \begin{equation} L = M_1 M_2\ldots M_{n-1} \end{equation} where $M_i = I + l_i e_i^{T}$ ...
3
votes
0answers
73 views

Efficient algorithm for a matrix porduct

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 ...
0
votes
0answers
35 views

Household reflectors to transform $T$ to upper triangular

Suppose you are given a nonsingular matrix $T\in\mathbb{R}^{n\times n}$. For example, if $n = 6$ then $$T = \begin{pmatrix} \alpha_1 & \beta_1 & 0 & 0 & 0 & 0\\ \gamma_2 & \...
0
votes
0answers
34 views

Compute an orthonormal bases for a full column rank matrix

Recall that given a full column rank matrix $A\in\mathbb{R}^{n\times k}$, we have discussed a reliable algorithm to compute an orthonormal basis of $\mathcal{R}(A)$ by computing the Household ...
3
votes
1answer
106 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 ...
0
votes
1answer
45 views

$LU$ Factorization of a nonsingular matrix with a particular pattern

Consider $S\in\mathbb{R}^{n\times n}$ whose nonzero elements have the following pattern for $n = 8$: $$\begin{pmatrix} 1 & 0 & 0 & 0 & \mu_1 & 0 & 0 & 0\\ 0 & 1 &...
1
vote
0answers
26 views

Schur complement of a matrix $A$

Let $A\in\mathbb{R}^{n\times n}$ and its inverse be partitioned $$A = \begin{pmatrix} A_{11} & A_{12}\\ A_{21} & A_{22}\\ \end{pmatrix},\:\: A^{-1} = \begin{pmatrix} \tilde{A_{22}} & \...
1
vote
1answer
46 views

Preconditioning of two step iteration for dense matrices

I would like to solve a dense linear system the form in python $$ L\left(\boldsymbol{x}\right):=\left[\gamma^+\left[\boldsymbol{A}+\frac{1}{2}\boldsymbol{B}^{-1}\right] +\gamma^-\left[\boldsymbol{A}-\...
0
votes
1answer
76 views

Parallel linear algebra without OpenMP

I have searched through the archives without success. Apparently, the question is simple: What linear algebra library can I use that is parallel (shared memory) but without OpenMP? As far as I've ...
1
vote
0answers
23 views

Optimal ordering in Jacobi SVD algorithm

In Jacobi SVD algorithm as given here every pair of columns of the matrix is orthogonalized until convergence. I want to know that how does the order of selection of the pair of columns affect the ...
6
votes
1answer
115 views

Is there an eigenvalue estimation method more accurate than Gershgorin's, which uses no multiplication?

Suppose I have a real symmetric matrix. I would like to tell wether it has at least $k$ strictly positive eigenvalues, but using only additions (no multiplications). Is there a method that I could use?...
0
votes
0answers
51 views

Pivoting in Crout algorithm

Consider the following code, found on wikipedia, that implements the Crout decomposition algorithm: ...
1
vote
1answer
83 views

Numerical method for solving a system with positive definite blocks

I have a system with below coefficient matrix $$ C = \begin{pmatrix} A & B^T \\ B & D \end{pmatrix},$$ where, $A$ and $D$ are square and positive definite. Furthermore, if $B$ be square, ...
3
votes
2answers
205 views

A numerical GMRES example

I'm having trouble understanding how GMRES works. I've read the part in Saad's book and a few others but still I am confused. Can someone provide me a numerical example to understand it better? Or if ...
5
votes
1answer
112 views

Solve for $C$ such that $C^{T}AC$ is banded of given width

Given a symmetric matrix $A$, the Lanczos algorithm outputs $C$ such that $C^{T}AC$ is tridiagonal. Is there a generalization of this such that $C^{T}AC$ is banded of specific width $w$? Note that $C$...
5
votes
0answers
189 views

Compute sparsity pattern of matrix product $A*A$

Suppose we have a sparse matrix $A$. Is there any faster way to compute just the sparsity pattern of $A^2=A*A$ (I do not actually need to know what exactly the nonzero value are) than to compute the ...
4
votes
0answers
109 views

Robust smoothers for geometric multigrid

I'm searching for robust smoothers for geometric multigrids. By robust I mean: Effective for high order approximations (say spectral element, spectral Discontinuous Galerkin), Parallel (suitable ...
1
vote
0answers
39 views

Fixing a near singular covariance matrix

Given a near singular covariance matrix, the standard method of 'fixing' it seems to be to add a small damping coefficient $c>0$ to the diagonal, which serves to bump all the eigenvalues up by this ...
0
votes
0answers
20 views

Algorithm to group Boolean functions by rotation similarity

I'm looking for a way to take the entire vector space of length n Boolean vectors and partition it into vectors that are the same up to a rotation of the entries. For example if n=3 the partitions ...
0
votes
2answers
20 views

Grouping Boolean vectors by similarity up to a rotation [duplicate]

I'm looking for an algorithm to take the entire vector space of length n Boolean vectors and partition it into vectors that are the same up to a rotation of the entries. For example if n=3 the ...
6
votes
2answers
171 views

Solving $A=B+AB$ without matrix inverse

I have a linear system of equations that can be expressed $$A=B+AB,$$ where $A$ and $B$ are real, symmetric matrices. I would like to solve for $A$ given $B$. At present, I solve for $A$ directly via $...
0
votes
1answer
57 views

Most efficient way to compute eigenvectors / values of this matrix?

I have a symmetric $ 3 \times 3 $ matrix $A$ and I need to compute the eigenvectors and eigenvalues of this. I know that I can use something like Lapack, but I also know that this can be computed ...
1
vote
1answer
50 views

Probability of reconstructing a word using c substrings from a random sample

Consider a voice recording split into it's phonemes as our sample $S=(s_1,...,s_k) \in \Omega = P^k$. The number of phonemes is $|P| = 40$. Then I have a word $w = (w_1,...,w_n) \in P^n$. I want to ...
1
vote
1answer
97 views

Numerical Solution of Matrix with Diagonal Elements of Highly Varying Order

I am trying to solve following set of equations: A(i,i-2)*u(i-2) + A(i,i-1)*u(i-1) + (A(i,i)+β(i) )*u(i) + A(i,i+1)*u(i+1) + A(i,i+2)*u(i+2)= B(i) + β(i) where i=1:1000000 If values of β varies(...
3
votes
0answers
43 views

Quasi Newton method for block diagonal Hessian

I am having an unconstrained optimization problem, where the Hessian is positive semidefinite and block diagonal. The function is strictly convex, hence, the curvature condition ($ s^{T}_{k}y_{k} > ...
1
vote
1answer
77 views

Efficiently rotate vector in 2D (and 3D)

I need to efficiently rotate a 2D (and 3D) vector in a CUDA kernel. I was thinking about generating random unitary rotation matrices. I don't need to know the angle, it just has to be randomly ...
5
votes
0answers
80 views

Galerkin FEM error when using even number of elements

Intro: I am developing a Galerkin FEM code in matlab - starting small with a simple 1D ODE. The equation I'm trying to solve is: $$a u_x = cos(x), x \in [0, 2\pi]$$ Which has a known exact solution $$...
0
votes
0answers
44 views

One parameter family of linear systems

Suppose $A\in\mathbb{R}^{n\times n}$ is dense and $D\in\mathbb{R}^{n\times n}$ is a diagonal and nonnegative. Given a fixed vector $b\in\mathbb{R}^n$, define the following function $f:\mathbb{R}\...
2
votes
0answers
26 views

SPECT reconstrction using MLEM

In Single-Photon Emission Computerized Tomography (SPECT) parallel beam reconstruction using Maximum-Likelihood Expectation–Maximization(MLEM), is it sufficient to scan the object around 180 degree? ...
0
votes
0answers
20 views

Ordering of eigenvectors to maximise trace of diagonalising matrix

I asked a similar question on the Mathematics stack exchange here without much success, so I thought I'd ask it with a more practical bent here. Suppose we have a Hermitian matrix $H$ with (for the ...
0
votes
0answers
55 views

Preconditioned Steepest Descent

For my program assignment I need to write a preconditioned steepest descent algorithm. I have psedo-code from my professor which is here: From reading in my text on this method they say that $P$ is ...
3
votes
1answer
63 views

Methods for solving rectangular, full-rank systems of equations — which is best?

Suppose I have a large, sparse, $m \times n$ matrix $A$, with $m \gt n$ and $\text{rank}(A) = n$. I wish to solve $Ax=b$. Suppose I know that $A$ has the following characteristics: $A$ is somewhat ...
0
votes
0answers
36 views

how to solve distance from overhead object to each of 3 ground-based radar receivers

How do I solve for the object distance to each receiver for three radar receivers on the ground, each the same distance from the other, and each receiving echoes, reflected from an object overhead, of ...
4
votes
1answer
64 views

Algorithm to decompose a sparse unitary matrix into a Kronecker product of smaller unitary matricies

Given some sparse unitary square matrix $A$ ($dim=2^n$ if it matters), is there an algorithm to decompose $A$ into a Kronecker/tensor product of smaller unitary matrices? In other words: decompose ...
4
votes
1answer
112 views

Power Iteration over Rayleigh Quotient Iteration?

It is a commonly known fact that the Rayleigh Quotient converges cubically (1), while the Power Iteration may converge slowly if the difference between the dominant and second-dominant eigenvalue is ...
8
votes
1answer
146 views

What is the state of the art algorithm for diagonalizing real symmetric matrices?

There are many methods for diagonalizing matrices, probably the most widely used is the combination of Householder transformations and the QR algorithm. Is there any superior method for diagonalizing ...
0
votes
0answers
116 views

Preconditioned Conjugate Gradient linear system solver in MATLAB

I have been trying to use the MATLAB's pcg() function to minimize an energy functional. Converting minimization problems to the solution of a linear system is ...
4
votes
1answer
103 views

Efficient way to generate a list of possible matrices (all integer components) with a determinant $V$

I have an interesting problem from my research that I have been struggling to solve. One part of the problem involves generating all possible matrices, where each set contains three integer vectors, ...