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

1
vote
0answers
23 views

Library for calculating determinants with Kronecker products

I need to calculate a determinant consisting of vectors, using the Kronecker product as product. As an example I would need to be able to calculate: $\left| \begin{array}{cc} ...
-2
votes
0answers
24 views

how could I find the minimum of computational cost of operation

Let the following statement: $$\phi=p_{xx}A^{-1}B+2p_x(A_x^{-1}B+A^{-1}B_x)+p(A_{xx}^{-1}B+2A_x^{-1}B_x+A^{-1}B_{xx}),$$ where $$A_x^{-1}=-A^{-1}A_xA^{-1}$$ and ...
2
votes
0answers
59 views

Using SVD to biorthogonalize left and right eigenvectors?

I have a set of left and right eigenvectors from an nonsymmetric eigenproblem, and I'd like to biorthogonalize them. I tried Gram-Schmidt, but this fails for most cases. I then read that the SVD is ...
1
vote
0answers
75 views

eigs routine in octave

I am using octave and observed a problem with the eigs-routine for non symmetric matrices. Using GNU octave version 3.8.1 the code below gives significant difference of eigenvalues although same ...
1
vote
0answers
58 views

Search Direction in Conjugate Gradient

Could you help me with a Conjugate Gradient question? In using CG to solve Ax=b, why is the search direction $p_{k+1}$ in CG chosen as a linear combination of the residual $r_k$ and previous direction ...
1
vote
1answer
91 views

Solve for a matrix given two vectors

I'm programming a beam finite element model by following a book (Nonlinear Finite Element Analysis of Solids and Structures Volume 2, in case you're wondering!). I've come across the following ...
2
votes
0answers
42 views

Algorithm code for Drazin and Bott-Duffin inverse (Matlab or C)

I could find the common Moore-Penrose algorithm, but I couldn't find the Drazin or the Bott-Duffin generalised inverse, except for some very specific cases, useless for my studying purposes. Is there ...
5
votes
0answers
52 views

Find the solution of linear equation using Wiedemann/ Krylov method

Let given $M =$ 1 0 1 0 1 1 1 1 1 and $b =$ 1 0 1 How to find the solution $x_3$ where ...
1
vote
1answer
59 views

Numerical eigenbasis for a unitary matrix

Do you know what numerical software computes an eigenvector basis for a unitary matrix? Say I have a unitary matrix $U$. If its eigenvalues are simple (no multiplicities), then for instance Matlab ...
2
votes
0answers
37 views

In-place QR update: deleting a column

Background I'm trying to do an update to a "thin" QR decomposition ($A = QR$, where $Q$ is $\mathbf{R}^{m,n}$, the first few columns (up to the matrix rank) of an orthogonal matrix and $R$ is upper ...
2
votes
2answers
187 views

Eigenvectors: MATLAB vs LAPACK DGGEV or DGGEVX

If we call LAPACK DGGEV or DGGEVX routines for two badly-conditioned matrices in a C++ code, will we get the same eigen-values ...
3
votes
1answer
130 views

Why is my MATLAB code for back-substitution slower than the backslash operator?

I wrote the code below to invert an upper triangular matrix, avoiding any possible multiplication/subtraction by zero. It just uses $\frac{1}{6}n^3+\ldots$ flops instead of $n^3+\ldots$ flops. ...
2
votes
3answers
125 views

Is it better to do normalization after all orthogonalization in Gram-Schmidt process?

In Gram-Schmidt process, is it better to do normalization after orthogonalization of all the vectors in a basis, or to normalize each new vector immediately after it is created, from computational ...
5
votes
1answer
89 views

Solving Generalization of Saddle point problem

I am interested in knowing if there is a generalization of the Uzawa iteration for the linear problems of the form $$\left[ \begin{array}{ccc}A& B^T&0\\ B&0&C^T\\ 0&C&0 ...
0
votes
1answer
79 views

Diagonalize a circulant-plus-rank-one matrix

Suppose $A=uv^T$ where $u$ and $v$ are non-zero column vectors in $\mathbb{R}^n, n\geq 3$. $\lambda=0$ is an eigenvalue of $A$ since $A$ is not of full rank. $\lambda=v^Tu$ is also an eigenvalue of ...
1
vote
1answer
100 views

writing linear system in matrix form

I have the following set of linear equations $$a_{m+1,n}+a_{m-1,n}+m(a_{m,n+1}+a_{m,n-1})+(m^{2}+n^{2})a_{m,n}=f_{m,n}$$ Here $m$ and $n$ run from 1 to $N$, so there are $N^2$ equations for the ...
3
votes
1answer
101 views

Get symmetric Finite Difference matrix in non Laplacian settings

I would like to solve a system of differential equations $u+\nabla(\nabla\cdot u)=f$ or in more detail $a+\partial_t^2a+\partial_t\partial_xb+\partial_t\partial_yc=f$ ...
5
votes
1answer
144 views

Least-squares for a diagonal matrix

This is a follow-up to a different question I asked with more detail. For $v\in\mathbb{R}^n$, denote $D_v\in\mathbb{R}^n$ as the diagonal matrix with elements in $v$. Given a "tall" matrix ...
3
votes
3answers
136 views

Large overdetermined system of linear equations

I'm looking for a method to solve a large overdetermined system of linear equations in a least squares sense. The matrix is dense. I'd like to use a method that works even with limited memory (we ...
4
votes
2answers
96 views

Solving “Hadamard systems”

Suppose we have two matrices $A$ and $B$ (we can assume they're symmetric; if absolutely necessary I think they may be positive definite). Then, is there any technique for solving $$(A\circ B)x=b,$$ ...
0
votes
2answers
88 views

How to determine the number of c points in algebraic multi grid

I am trying to write an algebraic multi-grid solver (in c++). At a given level I determine which nodes are c-points and which nodes are f-points (where the total number of c and f points equals the ...
2
votes
0answers
71 views

Does it matter if I use principal component analysis on the transpose instead of the original matrix?

My data set is a 60x10 matrix. I performed principal component analysis of this matrix with matlab using the princomp(AdjustedData) after I adjusting my original data set by subtracting the mean of ...
4
votes
2answers
104 views

Schur(QZ) Decomposition Differences

I am having issues understanding why different languages are producing different answers for the Schur(QZ) decomposition. I am working on writing some old stuff from Matlab into Julia and Python and ...
1
vote
0answers
28 views

linear solution of curve fitting on multiple linear functions differing by a multiplier

I am facing the following problem. I know nonlinear least squares can provide a solution but I am wondering if a linear way to solve this data fitting problem may exists. This is my input dataset: ...
1
vote
1answer
104 views

What is the algorithm that matlab used in its built-in function 'pca'?

Do anyone know what is the algorithm that MATLAB used in its built-in function "pca"? I have the following data set: 148.9820 55.8438 210.2150 149.3030 56.8891 208.4280 151.4400 ...
3
votes
1answer
169 views

Parallel Gram-Schmidt algorithms

I've heard that classical Gram-Schmidt is more amenable to parallelization than modified Gram-Schmidt; apparently the reason has something to do with "level 2 BLAS", which I'm not familiar with. ...
1
vote
1answer
129 views

Translate a 3D point along a heading

I need to translate a point (P1) in 3D a certain amount, call it stepSize, along a vector described by a heading composed of ...
0
votes
0answers
62 views

LU Factorization update when adding columns

I am looking for a way to update the LU factorization of a general $m \times n$ matrix after adding a column to the matrix. I have to iterate this procedure so I will begin with a matrix that is $m ...
0
votes
0answers
21 views

Updating the LU Factorization [duplicate]

I am looking for a way to update the LU factorization of a general m×n matrix after adding a column to the matrix. I have to iterate this procedure so I will begin with a matrix that is m×1 and go all ...
5
votes
2answers
186 views

Dimensionality reduction of the domain of f(x)

I'm wondering if there is something analogous to a PCA for data sets where there is a dependent variable. (Though I am interested in any method of dimensionality reduction, PCA is just an example.) ...
0
votes
0answers
49 views

Principal Components Analysis Not Behaving as I Expect it to

I have a bunch of points in $\mathbb{R}^3$ that I would like to translate and rotate so that their center is at the origin and the variance along the $x$ and $y$ axes are maximal (greedy, and in that ...
6
votes
1answer
141 views

Efficient RQ decomposition

I have an upper trapezoidal matrix stored in column major format. That is, my matrix looks like this: I'd like to RQ decompose it, and store Q in the rectangular part of my upper trapezoidal ...
2
votes
1answer
198 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
88 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
93 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
95 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
232 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
71 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
57 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
116 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
78 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
60 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
207 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
159 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
110 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
75 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
159 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, ...
15
votes
5answers
539 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
184 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
63 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 ...