Questions on the algorithmic/computational aspects of linear algebra, including the solution of linear systems, least squares problems, eigenproblems, and other such matters.
4
votes
2answers
118 views
Computing matrix exponential with PETSc/SLEPc
My question is basically the best way to compute the matrix exponential in parallel. I've got a sequential code written in Python (https://github.com/hadsed/AdiaQC) and I need to extend it to run in a ...
7
votes
1answer
55 views
Preconditioning and effects on precision of solution of LSE
In my courses on numerical analysis, I have been taught that the main and principal motivation for preconditioning linear systems of equations is to increase the convergence rate of iterative solvers ...
4
votes
1answer
94 views
generalized eigenvalue problem
I need to solve a real generalized eigenvalue problem
$Ax= \lambda Bx(*)$
A and B are calculated from equations below:
$$A=\sum_{i,j=1}^{N}W_{ij}(K_{i}-K_{j})\beta\beta^{T}(K_{i}-K_{j})^{T}$$
...
4
votes
1answer
61 views
$\mathbf{UDU}^\top$ decomposition routines in LAPACK/Eigen?
I would like to compute the decomposition of a real symmetric positive definite matrix $\mathbf{A} = \mathbf{UDU}^\top$.
LINPACK seems to have it as DSIFA, but I ...
2
votes
1answer
81 views
SVD and Lanczos method
I have been given the task of implementing SVD using the Lanczos method.
As I understand it only helps to find eigenvalues using a sequence of tridiagonal matrices, but i don't know how to apply this ...
3
votes
2answers
111 views
Quadratic Programming: Quadprog
Given a real-rectangular matrix $S$ and inorder to solve this simple quadratic programming problem:
Minimize $w'S'Sw = ||S w||^2$ over $w$ subject to
$e^Tw = 1$ and $w \geq 0$
using a solver I ...
7
votes
2answers
92 views
Is the maximum/minimum principle of the heat equation maintained by the Crank-Nicolson discretization?
I'm using the Crank-Nicolson finite difference scheme to solve a 1D heat equation. I'm wondering if the maximum/minimum principle of the heat equation (i.e. that the maximum/minimum occurs at the ...
1
vote
1answer
66 views
LU Decomposition with memory-mapped matrices
I have a ~4.12 Tb structured relatively-sparse matrix dataset (about 8% of the matrix entries are non-zero) that i want to apply an LU decomposition, however, given the size of it, loading it in ...
3
votes
1answer
66 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 = ...
2
votes
3answers
147 views
Can anyone recommend a library in C++ which has the most efficiency in doing sparse matrix operations under Windows system
I have download sparselib++, but it seems that it can't be complied in Windows, only in Unix, I don't know. So can anyone recommend some which can be used in Visual Studio 2010? I want to do a large ...
0
votes
1answer
63 views
Convergence of GMRES
From what I understand the GMRES method is (using Arnoldi Iterations/Modified Gram-Schmidt):
The first vector of the Krylov subspace span of A is the normalized vector $\frac{\vec b - A\vec x_0} {|| ...
2
votes
1answer
52 views
Solving the elliptic eigenproblem with periodic boundary conditions
Given an energy level $\mu$, I'm looking to calculate the eigenvectors corresponding to the time-independent Schrodinger operator on the torus (that is, periodic boundary conditions)-- $H = -h^2 ...
1
vote
1answer
66 views
Solving a “generalized” linear programming problem
I need to solve the follwing constrained optimization problem. Is there any MATLAB
toolbox which can solve this.
minimize $A_1 x$
subject to $A_2 x = 0$
$A_1$ is an $m\times n$ matrix where $m ...
1
vote
0answers
28 views
On solution of a class of discrete-time Lyapunov equation for systems with multiplicaitve noise
Let's consider the following equation
$$X=F_{1}XF_{1}^{T}+...+F_{p}XF_{p}^{T}+C$$
where $p$ is an positive integer and $C$ is a known positive semidefinite matrix. If we augment $F=[F_{1}...F_{p}]$ ...
0
votes
1answer
83 views
Issues with solving large sparse linear equations
I have some issues solving sparse linear equations Ax = b
My matrix A is sparse with dimension of 5 million by 5 million. Actually, it is a combination of two matrices. One is tridiagonal and the ...
10
votes
2answers
234 views
How is the SVD of a matrix computed in practice
How does MATLAB, for instance, calculate the SVD of a given matrix? I assume the answer probably involves computing the eigenvectors and eigenvalues of A*A'. If ...
3
votes
2answers
62 views
Parameter selection difficulty in GMRES
we are using a first order implicit finite volume code for simulation of incompressible flows. At its core, the code utilizes a (non-preconditioned) GMRES method for solving linear systems given in ...
1
vote
0answers
78 views
Block Cyclic Reduction implementation help, matlab
I'm sorry if this question is too open ended, or not appropriate for the forum. If so I'll remove the post and see if I can refine my problem, but at the moment I'm totally stumped.
I was hoping ...
1
vote
1answer
43 views
How does one implement the divide-and-conquer eigenvalue algorithm for a complex Hermitian matrix?
Let $T$ be a real symmetric tridiagonal matrix. Then the divide-and-conquer eigenvalue algorithm, as detailed in any standard text, goes by subdividing
$$ T = \begin{bmatrix} T_1 & 0 \\ 0 & ...
4
votes
1answer
82 views
Sparse matrix implementation of the Kalman Filter?
I have a Kalman Filter based modelling code that I have developed for a near-real time regional ionospheric mapping application. The code assimilates data from different sensors into a map (described ...
2
votes
2answers
112 views
Efficient computation of Markov chain transition probability matrix
Consider a continuous Markov chain $X=(X_t)$ on a finite state space and let $Q$ be the (given) transition rate matrix. This matrix is very sparse, with non-zero values on 3 diagonals only (so from ...
3
votes
1answer
42 views
explicitly forming coarse matrices with polynomial smoothing AMG
I've been reading about the algebraic multigrid algorithm and came across polynomial smoothers in this paper. It's my understanding that usually the coarse-level matrices
$A_H = I_h^HA_hI_H^h$
are ...
8
votes
2answers
119 views
Solving a linear system with matrix arguments
We're all familiar with the many computational methods to solve the standard linear system
$$
Ax=b.
$$However, I'm curious if there are any "standard" computational methods for solving a more ...
1
vote
1answer
60 views
Ground state eigenvector different for different eigen solvers (differs by negative sign in the elements). Does it matter?
Here is some code that hopefully clearly illustrates what I'm doing:
...
3
votes
2answers
138 views
Finding A and X such that AX = 0, X is positive non-zero, and A is sparse
I apologize if this is a naive question. I'm trying to create some boostrap data for a system of linear, ordinary differential equations at steady state.
Since the equations represent the ...
1
vote
1answer
76 views
4th order Padé scheme formula derivation
I am trying to derive the formula of the 4th order Padé scheme that passes through the points $x_i$, $x_{i-1}$ and $x_{i+1}$
$$\Big(\frac{\partial\phi}{\partial x} \Big)_i = ...
5
votes
3answers
136 views
computing the determinant of a dense nonsymmetric 100x100 matrix having very big and very small eigenvalues
The motivation for my question is the following: in one of Project Euler questions there is a need to count the spanning trees of a rectangular grid graph of dimension 100x500. By the Matrix-Tree ...
1
vote
0answers
27 views
Pixel-To-Angle Transformation in Camera Image
I'm trying to localize points I see in a camera image in terms of azimuth and elevation and match points between shots.
Individual shots should differ only in rotation around the camera's center ...
6
votes
1answer
45 views
Bad scaling versus collinearity
I was trying to solve a linear system:
$$
\mathbf{A}\mathbf{x} = \mathbf{y}
$$
but the conditioning number was quite bad (around $10^{17}$).
I thought that the system was singular, but after scaling ...
7
votes
1answer
123 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 ...
6
votes
3answers
128 views
Numerically stable explicit solution of small linear system
I have an inhomogeneous linear system
$$
Ax=b
$$
where $A$ is a real $n\times n$ matrix with $n\leq 4$. The nullspace of $A$ is guaranteed to be of zero dimension so the equation has a unique ...
3
votes
5answers
269 views
MATLAB's CVX Package to minimize $\mathrm{trace}(S)+\mathrm{trace}(S^{-2})$
I would like to minimize in matlab the function $f(S)=\mathrm{trace}(S)+\mathrm{trace}(S^{-2})$ where $S\in \mathcal{M}_{m,m}$ symmetric positive definite which is definitely convex function.
So I ...
3
votes
1answer
33 views
How to obtain the minimum set of variables required in a model to produce accurate estimation?
I have a system which I assume is linear. I have a matrix $A$ of which each row is a coefficient of a unknown variables in vector $x$. I have vector $B$ which contains the result of each $Ax$.
...
2
votes
0answers
28 views
Left and right eigenspaces of product of grammians
I solve the Lyapunov equations :
$$ A W_C E^T + E W_C A^T + B B^T = 0 $$
$$ A^T W_O E^T + E W_O A + C^T C = 0 $$
to obtain $ W_C $ and $W_O$. My aim is to get the left and right eigenspaces of $W_C ...
8
votes
0answers
69 views
Specialized methods for symmetric tridiagonal generalized eigenvalue problems
I have to solve generalized eigenvalue problems $Ax = \lambda Bx$ where $A$ and $B$ are both tridiagonal, $B$ is symmetric positive definite and real, but $A$ is only complex symmetric (not definite ...
3
votes
1answer
54 views
simplifying a product of a determinant and an inverse of a (nearly) singular matrix
Given two square matrices, $A$ and $B$, I need to calculate the product $tr(A^{-1}B)\times detA$. The catch is that $A$ is singular --- more precisely, it depends on some parameter $t$, such that it's ...
4
votes
1answer
62 views
Sparse LU for block-sparse matrices
I frequently need to solve linear systems with sparse matrices of moderate dimension (say a few thousand). These matrices are composed entirely of small dense blocks (typically 5-10 in dimension), and ...
4
votes
2answers
92 views
Complex least-squares problem
Having a matrix $\mathbf{A} \in \mathcal{C}^{m\times n}$ I solve following least-squares problem $$Re(\mathbf{A}^H \mathbf{A})x=Re(\mathbf{A}^H\mathbf{b}).$$ If the matrix $\mathbf{A}$ was a real ...
5
votes
2answers
224 views
Computing the pseudoinverse of a 3x3 matrix
I need to compute the (Moore-Penrose) pseudoinverse of fixed-size 3x3 matrices. I would prefer to have a simple method without bringing in the industrial strength machinery of Lapack. Are there any ...
8
votes
5answers
271 views
Repeatedly solving $A\mathbf{x} = \mathbf{b}$ with same $A$, different $\mathbf{b}$
I am using MATLAB to solve a problem that involves solving $A\mathbf{x}=\mathbf{b}$ at every timestep ($\mathbf{b}$ changes with time). Right now I am accomplishing this by MATLAB's ...
4
votes
3answers
123 views
implicit vs. explicit domain decomposition methods
I've been working on a finite element code on unstructured methods, which I've parallelized using the Schur complement method. Here's a summary of how I did it:
Assign each triangle of the mesh to a ...
2
votes
1answer
66 views
Algorithms for Compressed Sparse Rows
Is there a general survey on the basic algorithms for the compressed sparse rows format (like transposition, multiplication, addition, ...)? While it is not hard to write effective algorithms for ...
3
votes
1answer
65 views
How to pick a basis for the result of a non-linear function given a basis for its argument
I am trying to represent the result of a non-linear function in a small basis, given another small basis that does a good job a representing the argument of the function.
More specifically, there ...
3
votes
1answer
80 views
Numerical Methods for minimizing a Non-Differentiable Convex Function of Several Variables
I have a multi-variable convex continuous function which is not differentiable. I am interested to know about different numerical techniques, possibly also references to them, used for this.
Read ...
5
votes
1answer
127 views
Eigenvalue Decomposition of Hermitian Matrix in Scala
I'm working on helping my friend create code to perform the Overlap Dirac Operator and have come across one part that I'm not sure what to do.
I need to compute the Eigenvalues and corresponding ...
3
votes
1answer
108 views
restriction and interpolation in multigrid method
I need detailed explanation of the formula below
A2=I1*A1*I2
I suppose this formula computes matrix A2 on a coarse grid and here A1 is original matrix on fine ...
4
votes
1answer
108 views
Simple Lanczos algorithm code to obtain eigenvalues and eigenvectors of a symmetric matrix
I would like to write a simple program (in C) using Lanczos algorithm. I came across a Matlab example which helped me to understand a bit further the algorithm, however from this piece of code I can't ...
2
votes
1answer
59 views
Introduction for (numerical) linear algebra of random variables
I am in search of an introduction into numerical linear algebra - or, at least, pure linear algebra - that treats the case when the input data are random variables.
A typical application would be to ...
2
votes
1answer
80 views
cholesky factorization of block matrices
I have a block matrix (either 2x2 blocks or 3x3 blocks) which is the covariance matrix for a joint space of two or three multivariate normal variables. ie
...
6
votes
2answers
176 views
Finding the distribution (histogram) of eigenvalues for large sparse matrices
Are there any existing programs that are able to compute the (approximate) distribution of eigenvalues for very large (symmetric) sparse matrices?
Note that I do not need the eigenvalues themselves, ...
