A message from our CEO about the future of Stack Overflow and Stack Exchange. Read now.

# Tag Info

20

Your matrix is of size 15,000 x 15,000, so you have 225M elements in the matrix. This makes for roughly 2GB of memory. This is much more than the cache size of your processor, so it has to be loaded completely from main memory in every matrix multiplication, making for approximately 100GB of data transfers, plus what you need for the source and destination ...

17

I'm going to disagree with some of the other answers and say that I believe that figuring out how to use LAPACK is important in the field of scientific computing. However, there is a large learning curve to using LAPACK. This is because it is written at a very low level. The disadvantage of that is that it seems very cryptic, and not pleasant to the senses....

16

Matrix exponentials of skew-Hermitian matrices are cheap to compute: Suppose $A$ is your skew-Hermitian matrix, then $iA$ is Hermitian, and via zheevd and friends you can get the decomposition $$iA = U \Lambda U^H,$$ where $U$ is the unitary eigenvector matrix and $\Lambda$ is real and diagonal. Then, trivially, $$A = U (-i \Lambda) U^H.$$ Once you ...

15

There are a number of issues in your question. Do not use Gaussian Elimination (LU factorization) to calculate the numerical rank of a matrix. LU factorization is unreliable for this purpose in floating-point arithmetic. Instead, use a rank-revealing QR decomposition (such as xGEQPX or xGEPQY in LAPACK, where x is C, D, S, or Z, though those routines are ...

15

As far as I know, Lapack is the only publicly available implementation of a number of algorithms (nonsymmetric dense eigensolver, pseudo-quadratic time symmetric eigensolver, fast Jacobi SVD). Most libraries that don't rely on BLAS+Lapack tend to support very primitive operations like matrix multiplication, LU factorization, and QR decomposition. Lapack ...

13

When you use ZGELSS to sovle this problem, you're using the truncated singular value decomposition to regularize this extremely ill-conditioned problem. it's important to understand that this library routine is not attempting to find a least squares solution to $Ax=b$, but rather it is attempting to balance finding a solution that minimizes $\| x \|$ ...

10

PETSc uses BLAS for a few vector primitives, but these are generally limited by memory bandwidth and there isn't much variance in "optimization", so it tends not to make much performance difference. It also uses Lapack for some analysis such as Lanczos or Arnoldi estimates of eigenvalues and singular values, but these are generally not performance-sensitive....

10

Trivial answer for square $A$: use dgesvx which solves also for $A^T x = b$ when TRANS = 'T'. Please note that with BLAS or LAPACK you hardly have to transpose (swapping elements in memory) a matrix: most of the subroutines have a TRANS argument to accommodate for operation on the transpose matrix or on a matrix stored with a different memory layout. (...

9

I can't answer the second half of your question as far as other implementations out there but I can provide some insight as to the challenges. For reference, I personally used ViennaCL on a nVidia GTX 560 Ti with 2GB of memory for my benchmarks. Over serial code on a mid-range i5, I saw speed-ups for dense matrix multiplications of approximately 40x. For ...

9

Jed Brown has already pointed this out in the comments to the question, but there is really not very much you can do in usual double precision if your condition number is large: in most cases, you will likely not get a single digit of accuracy in your solution and, worse, you can't even tell because you can't accurately evaluate the residual corresponding to ...

9

OpenBlas is quite fast, so you can link it to LAPACK. Have you tried precompiled version of LAPACK/BLAS from your CPU vendor? For example AMD ACML (free) or Intel MKL (free on linux for non-commercial and non-academic use)? You simply need to unpack and run install file. In my opinion the only advantage of using ATLAS is then when you use some unusual CPU. ...

9

There's no reason to append a row of 1's. You should just perform a rank-revealing QR factorization (like with routine SGEQP3) on $A^T$, and the last column of $Q$ should be in the nullspace. This has the added advantage that the relative magnitude of the last element on the diagonal of $R$ gives you some idea of how singular the solution is. Even better ...

9

This can be interpreted as summing over an index of a tensor when the vector $x$ is reshaped into a box of numbers instead of a list. In particular, if $X$ is the $d\text{-by-}d$ folded version of $x$, then the operation you are doing is, \begin{align} Dx &= \mathrm{vec}\left((I \otimes \mathbf{1})\mathrm{vec}(X)\right) \\ &= \mathrm{vec}(\mathbf{1}^...

8

Short answer: nothing more than $U_{ii} = 0$, i.e. that your computed $U$ factorization is exactly singular. xGETRF is not safe as a rank revealing factorization, so I would not draw any conclusion, apart from the fact that $A$ is ill-conditioned and no solution to $Ax=b$ can be safely computed. Information on rank and null space should be derived (if ...

8

I would first think really hard about whether or not the matrix is really completely arbitrary: Is there any transformation that would make it Hermitian? Does the physics guarantee that the matrix should be diagonalizable (with a reasonably conditioned eigenvector matrix)? If it turns out that there really isn't any symmetry to exploit, then you should ...

8

Due to the wording of your question, I am assuming that your matrix is square. LAPACK's SVD routines, such as zgesvd, essentially proceed in three stages for square matrices: Implicitly computing unitary matrices $U_A$ and $V_A$, as products of Householder transforms, such that the general matrix $A$ is reduced to a real, upper bidiagonal matrix $B := U_A^H ... 8 There is a neat trick I have recently learned from this paper. You start doing rank-revealing QR, and stop after the first$k$Householder reflections, when you have a matrix of the form $$\begin{bmatrix} R_1 & R_{12}\\ 0 & R_{22} \end{bmatrix},$$ with$R_1$triangular of size$k\times k$, and$R_{22}$typically not triangular (since we stopped ... 7 To build on what Jack has said, the standard approach that seems to be used in software (like EXPOKIT, mentioned in your earlier question) is scaling-and-squaring followed by Padé approximation (Methods 2 and 3) or Krylov subspace methods (Method 20). In particular, if you're looking at exponential integrators, you'll want to consider the Krylov subspace ... 7 There are two issues at hand here: Is$A$dense or sparse? Do you have the same software stack as MATLAB's internal libraries? Dense or sparse? MATLAB no longer explicitly mentions the LAPACK routines it calls to obtain a QR factorization if$A$is dense. If the information in the documentation for MATLAB R2008b also holds for later releases, then MATLAB ... 7 numpy.linalg.svd is a wrapper around {Z,D}GESDD from LAPACK. LAPACK, in turn, is very carefully written by some of the world's foremost experts in numerical linear algebra. Indeed, it'd be very surprising if someone not intimately familiar with the field would succeed in beating LAPACK (either in speed or accuracy). As for why QR is better than Gaussian ... 7 The simplest/fastest way to solve ill-conditioned problems is to increase precision of computations (by brute force). Another (yet not always possible) way is to re-formulate your problem. You might need to use quadruple precision (34 decimal digits). Even though 20 digits will be lost in a course (because of condition number) you will still get 14 correct ... 7 The repository package is not safe to use with threading due to the way it was compiled. I reported the bug on the Lapack forum, but it will take a long time for workarounds or solutions to trickle down into the repository. If you compile it yourself, be sure to add the "-frecursive" to gfortran. 7 Took me quite a while to figure this out and as usual it becomes obvious after you find the culprit. After checking the problematic cases reported in David S. Watkins. A case where balancing is harmful. Electron. Trans. Numer. Anal, 23:1–4, 2006 and also the discussion here (both being cited in arXiv:1401.5766v1), it turns out that matlab uses the ... 7 I usually resist telling people what I think they should do rather than answering their question but in this case I'm going to make an exception. Lapack is written in FORTRAN and the API is very FORTRAN-like. There is a C API to Lapack that makes the interface slightly less painful but it will never be a pleasant experience to use Lapack from C++. ... 7 Here's another answer in the same vein as the above. You should look into the Armadillo C++ linear algebra library. Pros: The function syntax is high-level (similar to that of MATLAB). So no DGESV mumbo-jumbo, just X = solve( A, B ) (although there is a reason behind those oddly-looking LAPACK function names...). Implements various matrix decompositions ... 7 The problem, of course, is that computing the true rank (e.g., via a QR decomposition) is not really any cheaper than computing a low-rank representation of the matrix. The best you can probably do is to use a randomized algorithm to find low-rank approximations. These can, at least in theory, be significantly faster than working on the entire matrix ... 7 It's the blocked variant of Householder-QR that is driving this design. If you look in Golub and Van Loan's book (Ch 5.2 or so) they talk about how k-iterations of the algorithm can be blocked together by accumulating the individual reflectors into a rank-k reflector of the form$\mathbf I + \mathbf W \mathbf Y^{\mathrm T}$, where both$\mathbf W$and$\...

6

See Leslie Foster's page on rank-revealing software. See also this LAPACK Working Note analyzing failures of rank-revealing QR xGEQP3. You should be able to find out what routines MATLAB uses by setting breakpoints in a debugger and examining the stack. Last I looked, admittedly several years ago, MATLAB used shared libraries, in which case the symbol names ...

6

You are using a reference implementation that does partial pivoting. Tridiagonal solves do very little work and do not call into the BLAS. It is likely slower than your code because it does partial pivoting. The source code for dgtsv is straightforward. If you will solve with the same matrix multiple times, you may want to store the factors by using dgttrf ...

6

In my experience, the best way to use blas/lapack on recent versions of ubuntu is to use the packaged openblas. For what it's worth, I mostly use blas/lapack through python numpy/scipy, and using openblas speeds up some of the linear algebra by like 200x vs. the default. I've tried using custom ATLAS, but it was a huge pain and didn't give much if any ...

Only top voted, non community-wiki answers of a minimum length are eligible