Hot answers tagged

65

Except for code which does a significant number of floating-point operations on data that are held in cache, most floating-point intensive code is performance limited by memory bandwidth and cache capacity rather than by flops. $v$ and the products $Av$ and $Bv$ are all vectors of length 2000 (16K bytes in double precision), which will easily fit into a ...


21

Your code is limited by memory bandwidth. For trivial math, it's often better to count memory accesses rather than flops. You'll get the following table: operation memory reads/writes matrix + matrix 3n² matrix * vector 2n²+n (if vector is not cached) matrix * vector n²+2n (if vector is only read once) vector + vector ...


12

If you can compute products with $A$ and $A^T$, as you specify in a comment, you can run the classical sparse SVD algorithms such as scipy.sparse.linalg.svds, Matlab's svds, or Julia's Arpack.svds, which are based on Lanczos bidiagonalization. They are designed to compute singular values, and are likely to be more robust than a minimization routine coded by ...


7

You are effectively asking how to compute $y=(\log M )^{-1}b$, where $M=e^A$ is the given matrix. There are several methods for computing $f(M)b$ without forming $f(M)$, and they are reviewed here. One general method is to use Cauchy’s theorem, $$y=\dfrac{1}{2\pi i}\int_\Gamma f(z)(zI - M)^{-1}b\,dz,$$ with $f(x) = 1/\log(x)$. $\Gamma$ is a contour that ...


7

I believe you can accomplish what you want efficiently using the recursive LU algorithm. In brief, recursive LU on a $M \times N$ matrix $A$ proceeds by partitioning the matrix into 4 blocks: \begin{align} \pmatrix{A_{11} & A_{12} \\ A_{21} & A_{22}} &= \pmatrix{L_{11} & 0 \\ L_{21} & L_{22}} \pmatrix{U_{11} & U_{12} \\ 0 & U_{22}}...


7

They are commonly called block matrices. You can create them with hstack, vstack, and block.


7

One part of the fundamental theorem of linear algebra is that the kernel/nullspace of $\mathbf A$ is orthogonal to the range of $\mathbf A^T$. By applying the $\mathbf Q \mathbf R$ decomposition to $\mathbf A^T$, you can generate the orthogonal projector $\mathbf P = \mathbf I - \mathbf Q \mathbf Q^T$. The vector $\mathbf P \mathbf x$ is what you're looking ...


7

In general, I agree with Chris's comment that using a compiled language with the allocation of the matrices on the stack can help significantly. Several possibilities if we are limited to Python and numpy: consider np.array vs np.matrix, it might happen that np.matrix is faster than np.array matrix-matrix product (it is unclear what you are using now, and ...


5

Orthogonal matrices are about as well-conditioned as you can get, but numerical errors still occur. One common error is loss of orthogonality. A fix for this could be to re-orthogonalize your columns after some number of multiplications. You can do this by just taking the QR decomposition of your matrix after some number of products and taking the orthogonal ...


4

I think you already know the answer- find the inverse of each of the small diagonal blocks. The details of how you might break the matrix up into the block in Python are something out of scope for this stackexchange group. As others have suggested you most likely don't really need the inverse of this matrix and could probably do everything you need with a ...


4

I don't think there is a way to display the used LAPACK function names natively during runtime using the interfaces provided by scipy.linalg. Depending on your goals you can: read the source code and deduce the logic from there. Unfortunately, this is not a runtime-use scenario, but a human analysis. fork your own version of scipy, add custom outputs (to ...


4

You might want to use the fact that: $$ ||A||_2=\sigma_\max(A) $$ where $\sigma_\max$ is the largest singular value. If you are interested in details, this Math SO question should be interesting. Thus, $$ ||A^{-1}||_2=\frac{1}{\sigma_\min(A)} $$ where $\sigma_\min$ is the smallest singular value. You certainly want to avoid the actual calculation of the ...


4

I have not quite understood what the dimensions of your matrices are (see comments). but the following code may give you a start. It is not efficient or particularly elegant, but you can improve from there: import numpy as np A = np.matrix('1 2; 3 4') B = np.matrix('1 2; 3 4') C = np.matrix('1 2; 3 4') result = (np.matmul(A.transpose(),np....


3

N * (G' .* N') A few releases ago Matlab introduced singleton expansion: in expressions like the one in parentheses the $n\times 1$ matrix G' is "upgraded" to an $n\times m$ matrix (with all columns equal) before the elementwise product. This new feature gives a cleaner way to implement some tricks that previously required bsxfun.


3

Test your code. "I don't know if my linear algebra routines work or not" is a problem you can solve easily. It is easy to check if you have computed the correct eigenvalues or not; just check that $VDV^{-1}=H$, or if you don't fancy the inversion $HV=VD$. If you are not sure that your individual pieces work, you are walking in the dark. Possibly a hot take: ...


3

I wanted to follow up with this question because I'm absolutely astounded by the performance improvement I was able to achieve using a C extension to python. I wrote a simple function in C that takes my (2,2,n) numpy array and performs the repeated matrix multiplication on it. I followed KobeGote's recommendation to hard-code the 2x2 matrix multiplication ...


3

Bareiss' algorithm is a better choice if you have to compute exactly determinants of integer matrices, as noted in the comments. When it comes to real (floating-point) or complex matrices, the main point to understand is that computing determinants is not a common task. In scientific computing, determinants are often a wrong choice: whatever you wish to do, ...


3

As I mentioned in my comment, due to that you are searching for a method based on Python and ideally available in NumPy or SciPy, my suggestion is to use scipy.linalg.expm. It's really easy to use basically you have the $X$ matrix and you just pass it to the scipy.linalg.expm class and it would give you the exponential of $X$ (i.e. $e^{X}$). Due to that you ...


3

As others have mentioned in comments, solving $A^TAx=A^Tb$ (the so-called least-norm solution) can be performed without explicitly forming $A^TA$ with the QR decomposition. It works like this: $$A=QR$$ where $Q$ has the interesting properties that $Q^TQ=I$ and $R$ is square and triangular. Using this information, the least squares system $A^TAx=A^Tb$ is ...


3

LAPACK has an implementation of the svd of a 2x2 triangular matrix. It appears to be very robust. The routine is XLASV2. To apply to a regular 2x2 matrix, you can simply apply a single givens rotation from the left/right.


3

The Frobenius norm of a binary matrix is the square root of the number of non-zero elements. Let the point (0.,...0) be the origin, and let's say the vec'd binary matrix elements are the coordinates of a point, with coordinates being zero or one. Then the Frobenius norm is the Euclidean distance from the origin to that point.


3

If your advection problem had Dirichlet of Neumann boundary conditions, the linear system would be tridiagonal and you could apply the Thomas algorithm. With periodic boundary conditions, however, we lose this. If c(x) is a constant independent of x, the matrix would be circulant and linear systems could be solved efficiently using FFTs. An even better ...


3

Another nice way to do this sort of thing is numpy.einsum, which allows you to express the multiplication in index notation: >>> A=np.array([[1,2,3],[4,5,6]]) >>> B=np.array([[1,2],[3,4]]) >>> C=np.array([[1,2,3],[4,5,6]]) >>> A.T@B@C array([[ 85, 116, 147], [113, 154, 195], [141, 192, 243]]) >>> np....


2

This is a binary knapsack problem which is known to be NP-hard. No efficient algo exists yet, but there are algos that can solve problem up to size of 400 variables (according to a paper published in 1999).


2

You are looking for Chamfer-Matching: it will yield a zero value at all center points where the pattern occurs. Ordinary template matching (aka cross correlation) does not work in your case, because that would also consider the white points of the template, not only the black points. Chamfer-Matching, instead, only sums the distance transform values of A at ...


2

I am expanding my comment into an answer. I don't think that the method is efficient, but I think that it can be used to obtain the matrix $A$ from $e^{A}$. We know that $$\frac{d e^{tA}}{dt} = e^{tA} A\, ,$$ so, we could use $$\left.\frac{d e^{tA}}{dt}\right|_{t=0} = A\, ,$$ if we can approximate the derivative $$\frac{d e^{tA}}{dt} \approx D(A)\, .$$ ...


2

Solution slightly modified from here W2 = N * bsxfun(@times, N, G).'; This works for N of size m x n, and G of size 1 x n.


2

I also suggest looking into the condition number estimators, which will (with some degree of [un]reliability) predict how effectively numerically singular the matrix is. In particular, "Spectral Condition-Number Estimation of Large Sparse Matrices" based on LSQR seems like an interesting choice. attracted my attention one day. I would also ...


2

@Federico Poloni 's fine answer states the impossibility of getting an exact yes/no answer using IEEE arithmetic. However, using interval arithmetic with outward rounding, it is possible to get a "not singular/don't know" answer. In particular, it may be possible to definitively conclude that the smallest singular value is strictly greater than ...


1

Here's my 2 cents. I would set up the following minimization problem $$ \pi(x) = \frac{1}{2} (Ax)^T(Ax) $$ If $A$ has eigenvalues which are zero, there will exist a nonzero $x$ such that $\pi(x)=0$. So, I would try computing the gradient of $\pi$ wrt $x$ and use a gradient descent algorithm to drive $\pi$ towards zero. If you get reasonably close to zero ($\...


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