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

2
votes
1answer
92 views

Comparing Eigenvectors, Mathematica vs. Matlab

I am trying to create the same out puts in Mathmatica and Matlab, however I am running into trouble aligning the eigenvectors with the eigenvalues, I think the Matlab is doing something slighly more ...
1
vote
3answers
65 views

rank-deficient NNLS

I want to find the minimum-norm solution to a rank-deficient least-squares problem, subject to positivity constraints, e.g. $$\min_x\ \|x\|^2 \quad s.t.\quad Ax = b,\ x \geq 0$$ where $A$ is large, ...
3
votes
1answer
83 views

How does matrix-matrix product scale with multiple CPUs?

These days, one can have 64 cores in a single node. I wonder how well the dense matrix-matrix product (SGEMM and DGEMM) scales ...
1
vote
1answer
41 views

Symbolic Computations with Block Matrices in Maple

Due to the derivation of an algorithm I tried to use Maple (9.5) to calculate some block matrix expression. Unfortunately Maple seems to ignore the assumption I set on the variable. Lets consider the ...
5
votes
2answers
210 views

Recommendations for symmetric preconditioner

Given symmetric, positive-definite matrix $A$ and its preconditioner $M^{-1}$. What's kind of preconditioner $M^{-1}$ which preserve the symmetry of $M^{-1}A$?
3
votes
1answer
54 views

State-of-the-art for active set optimization algorithms?

Given a problem like this: $$ \text{min } ||Ex-f|| \text{ s.t.}$$ $$ Gx \ge 0$$ $$ Cx = d $$ And assuming that the matrices are medium sized (dimensions in the low thousands) and dense, what's the ...
0
votes
0answers
34 views

Lapack : dgeev strange behaviour

I have a problem that is quite strange. I use c++ and dgeev Lapack to diagonalize and find the eigenvectors of a 36x36 real non symmetrical matrix. The documentation explains the following : if ...
5
votes
2answers
74 views

Sparse smallest eigenvalue problem on a linear subspace?

I am interested in methods for solving the optimization problem $$ \begin{array}{rl} \arg\min_x & x^T A x \\ \mathrm{s.t.} & x^T x = 1 \\ & Bx = 0 \end{array} $$ where $A$ is symmetric ...
7
votes
1answer
60 views

Do C++ matrix libraries translate compound vector operations to single loops?

I am trying to replicate Fortran90 array syntax using a C++ library. The libraries themselves are discussed at length in this question. They can all do something like this: ...
1
vote
1answer
53 views

Conditioning of matrix factorizations and square roots

For my application, I need factors $\tilde C$ so that $$ \tilde C{}^T \tilde C = C^TMC $$ where $C$ is a long matrix, i.e. $C$ has much more columns than rows, and $M$ is a small symmetric positive ...
5
votes
1answer
146 views

Testing if a matrix is positive semi-definite

I have a list ${\cal L}$ of symmetric matrices that I need to check for positive semi-definiteness (i.e their eigenvalues are non-negative.) The comment above implies that one could do it by ...
7
votes
1answer
118 views

Algorithm to calculate the exponential of an Hessenberg matrix

I am interested in computing the solution of a lage system of ODEs using a krylov method as in [1]. Such method involve functions related to the exponential (the so-called $\varphi$-functions). It ...
0
votes
1answer
64 views

solving underdetermined system of equations with a sparse matrix as input

I am using Matlab to solve Ax=b and my A is very large, sparse, binary and also rectangular. I saw the Matlab backlash \ operator help and it states that if A is rectangular then it will use the QR ...
3
votes
1answer
72 views

Threaded QR with column pivoting

My program needs to perform pivoted QR decomposition on tall (e.g. 1e9 by 100) matrices. I run into the bottleneck that the major computational time of my program is spent on doing serial pivoted-QR ...
4
votes
2answers
93 views

Spectral decomposition of symmetric matrix

What is a good direct method to compute the spectral decomposition / Schur decomposition / singular decomposition of a symmetric matrix? "Direct" means as in LU decomposition, Cholesky decomposition, ...
11
votes
3answers
376 views

20% performance penalty for a nice software design

I'm writing a small library for sparse matrix computations as a way to teach myself to make the best use of object-oriented programming. I've worked really hard on having a nice object model, where ...
0
votes
1answer
69 views

Flop counts for LAPACK symmetric eigenvalue routines DSYEV, DSYEVD, DSYEVX and DSYEVR

LAPACK has following 4 routines for calculating eigenvalues of a real symmetric matrix; namely DSYEV, DSYEVD, DSYEVX and DSYEVR (DSYEVR being the recommended one). If I were to calculate both ...
1
vote
2answers
59 views

Performance metrics to compare initial-boundary value problem solutions

I am comparing the performance several finite difference methods of solving an initial-boundary value problem. There are several dimensions to this comparison: Number of cells Number of timesteps ...
2
votes
3answers
228 views

Eigenvalues of a sparse banded nonsymmetric matrix from an elliptic operator

I have a sparse matrix coming from the discretization of a 3D elliptic PDE. The matrix is banded with seven non-zeros diagonals. The sparsity pattern of the matrix looks like this (the actual matrix ...
3
votes
3answers
282 views

Fast vector - “diagonal” matrix multiplication

Let $\mathbf{1}\in\mathbb{R}^d$ be a vector with all elements equal to $1$. Define: $$\mathbf{D} = \mathrm{diag}(\mathbf{1}^\top,\mathbf{1}^\top,\ldots,\mathbf{1}^\top) = \begin{bmatrix} 1 \cdots 1 ...
3
votes
2answers
154 views

Solving linear system with 6 equations and 22 unknowns for six of the unknowns

I am trying to find the solution for the M variables in the following system. \begin{equation} 0 = C_{b} M^{b}_{x} - M^{a}_{x} k_{2a} + M^{a}_{y} \left(\omega - \omega_{a}\right)\\ 0 = C_{a} ...
10
votes
2answers
211 views

preconditioning a krylov method with another krylov method

In method like gmres or bicgstab it could be attractive to use another krylov method as a preconditioner. After all they are easy to implement in a matrix-free way and in a parallel environment. For ...
5
votes
1answer
145 views

Matrix completion algorithm

I am trying to implement the algorithm presented in this paper which tries to recover a matrix that represent a less noisier dataset of the intensities of the pixels of a set of images. In this case ...
3
votes
1answer
85 views

Decompositions of sparse symmetric matrices and methods for solving large linear equations

I asked the same question on mo.se and it was suggested that scicomp would be a better forum for it. So here it is: I am writing code for solving linear equations of the form $$A_{n\times n}\cdot ...
2
votes
3answers
72 views

importing PARI libraries ( in C++) with Python

I noticed that Python does not have a good datatype for rational numbers, certainly not for algebraic numbers like $\tfrac{1 + \sqrt{2}}{3}$ or the real root of $x^3 - 5x + 7$. They have the ...
0
votes
0answers
45 views

3D surface reconstruction

I've been studying about Photometric Stereo to reconstruct a 3D surface from images but I'm kinda lost in the part of the Depth Map. I understand that the normal vector should be orthogonal to the ...
1
vote
0answers
44 views

Implementing the transition matrix for page rank

I'm trying to implement PageRank. I'm reading the description here: http://nlp.stanford.edu/IR-book/html/htmledition/markov-chains-1.html Everything is very clear to me, however I'm concerned about ...
1
vote
0answers
44 views

An efficient way to convert between MKL and Armadillo types

We use Armadillo a lot in our code but there are places where we prefer to use MKL directly from Vector Math Library. We have cx_vec i.e. vector of doubles. But the exp function of armadillo uses ...
0
votes
0answers
21 views

convergence of self-consistent solution- vector of large no. of points

I have a matrix equation $\left(\begin{array}{cc} \frac{\delta^{2}}{\delta x^{2}}+\mu & D(x)\\ D(x) & \frac{\delta^{2}}{\delta x^{2}}+\mu \end{array}\right)\left(\begin{array}{c} u(x)\\ v(x) ...
8
votes
2answers
177 views

Representing Eisenstein numbers without floats

I have a project where I need to use quadratic fields Specifically numbers of the form $a + b \sqrt{-3}$ with $a,b \in \mathbb{Q}$. For example here are the prime numbers in Eisenstein integers: ...
0
votes
1answer
53 views

Least squares fitting

I have the following equation I came across which was solved using least squares $x = \sum_{n=1}^{N} A_{n} y_{n}$ Where $x$ is a $m \times p$ matrix and $y$ would be of size $m \times p$ as well ...
2
votes
0answers
70 views

Inverse problem with a rank-1 update

I hope you can help me out with this. I have to find the solution x to an inverse system $$ x=A^{-1}b $$ This inverse problem is basically a least square problem with a rank-1 update. $$ ...
1
vote
1answer
77 views

LSA, SVD and the Frobenius norm

In Latent Semantic Analysis one uses the SVD to perform a dimensional reduction of the term-document matrix, via the Eckart-Young theorem. Now, the rank $k$ approximation obtained by E-Y is proven to ...
3
votes
3answers
121 views

Efficient solver for a symmetric tridiagonal system where the upper/lower diagonals are offset

I'm looking for an efficient way to solve a symmetric tridiagonal system $Mx = d$, where the upper and lower diagonals of $M$ are offset from the main diagonal by $k$ rows/columns: $$ \begin{bmatrix} ...
1
vote
0answers
53 views

Lanczos algorithm with thick restart on a dynamic matrix

According to a recommendation, this is a re-post of that. currently, I'm working on a way to compute the 2 biggest eigenvalues of a real, symmetric, huge and sparse matrix that changes a few entries ...
0
votes
0answers
110 views

How to solve singular non symmetric poisson equation with Neumann boundary condtions?

I am trying to solve 2D Poisson equations with Neumann boundary conditions. When the mesh is uniform, Poisson equation is singular and symmetric, so the method listed in Null Space Projection for ...
5
votes
2answers
149 views

Applying matrix square root inverse in matrix-free regime

Let $A$ be a large symmetric positive definite matrix, and suppose that we can efficiently apply $A$ and have a fast solver to apply $A^{-1}$, but we do not have access to the matrix entries for ...
3
votes
2answers
180 views

Methods for solving linear systems

This is such a basic topic but there are so many different methods proposed for solving a linear system of equations. I recently found a very good source but couldn't really make sense of all the ...
3
votes
3answers
179 views

How to solve a Linear Matrix Equation: AX-XA=B efficiently?

recently I have been working on solving some math problems using Fortran. There occurs to me that a linear matrix equation: $$ AX-XA=B $$ where $A$ and $B$ are known $n\times n$ matrices and $X$ is ...
2
votes
2answers
203 views

Good Finite Element Library for a small project

I'm currently working on this project and I have a basic structural analyzer that uses the finite element method. Essentially, I turn each block into a set of trusses, construct a stiffness matrix ...
5
votes
1answer
67 views

Computing eigendecomposition of a Hermitian matrix that is almost unitary

I have a dense Hermitian matrix that is approximately unitary, so it has eigenvalues that are $\sim \pm1$. I would like to compute all the eigenvectors corresponding to the $+1$ eigenvalue (not ...
5
votes
0answers
69 views

Wanted: sequences of linear systems for recycling Krylov solver analysis

In the solution of sequences of linear systems $$A_ix_i=b_i\quad\text{for}\quad i=1,2,\dots$$ with Krylov subspace methods, data can be recycled from already solved linear systems in order to speed up ...
2
votes
1answer
124 views

Weighted Gauss-Seidel Algorithm

In Jacobi method's Wikipedia article there's a section that describes Weighted Jacobi method: http://en.wikipedia.org/wiki/Jacobi_method#Weighted_Jacobi_method. I need to implement the Weighted ...
1
vote
2answers
109 views

What sparse linear programming solver it is better to use?

I have the following LP problem: $$ \min \limits_{\varepsilon, x_{1}, \ldots, x_{n}}f(\varepsilon, x_{1}, \ldots, x_{n}) = \varepsilon \;\;\;\;\; s.t. \;\;C x \geq 0, \;\; x_{i}^{0} - \varepsilon ...
2
votes
3answers
103 views

Iterating through a 3D triangle

Given a triangular plane formed by three points in R3 space {p1, p2, p3}, I want to iterate through all points on the triangle plane by using two variables, x0,y0, something like in this example: ...
1
vote
0answers
41 views

Sparse matrix factorization of a rank deficient matrix by decomposition into linearly independent components

I've got a little conjecture I need to prove for a theoretical result related to causal Bayes net search with latent variables under sparsity constraints. If you're interested in the application ...
3
votes
0answers
38 views

BLACS context value and multiple MPI communicators

I am reposting here a question previously asked in stackoverflow: as suggested in a reply to my question, scicomp could be a better place to obtain some useful comments/suggestions for my problem. I ...
1
vote
3answers
126 views

How to find the smallest positive eigenvalue of a large general system if they are all in +/- pairs of real eigenvalues

I used MKL Lapack SBGVX because it can solve for "selected" eigenvalues/modes (both positive and negative), thinking it would be efficient, but it is extremely slow when compared with Bathe's subspace ...
3
votes
1answer
100 views

spectral decomposition in Numpy, sign difference

I am trying to follow along with an example from a book, but I get seemingly different answers depending on which spectral decomposition function I use in Numpy. I am trying to transform the Matrix G, ...
4
votes
1answer
60 views

Does the covariance matrix in Least Squares depend upon the input data?

I had always assumed that the covariance matrix depends upon the amount and quality of your input data, but I am finding out that this is not the case. Is this true? We want to fit $f(t) = ...