Hot answers tagged

24 votes
Accepted

Complexity of matrix inversion in numpy

(This is getting too long for comments...) I'll assume you actually need to compute an inverse in your algorithm.1 First, it is important to note that these alternative algorithms are not actually ...
9 votes
Accepted

Alternatives to numpy.einsum

Pyhton 3.5 will introduce a new operator @, which was proposed by NumPy devs to be the matrix multiplication operator. You may want to read PEP 465, but the syntax ...
  • 255
9 votes
Accepted

Poor SVD reconstruction of singular matrix

Algorithms for the SVD, as more or less every classical linear algebra algorithm based on orthogonal transformations, are normwise backward stable, i.e., it should be guaranteed that $\frac{\|USV^* - ...
8 votes

Why `import numpy as np` for scientific computing?

I would say, that it can be explained by the following famous programming principle: Explicit is better than implicit Usually, that is applied to types; however, it can be applied to namespaces (...
  • 8,521
8 votes
Accepted

Fastest Way to Mutiply $10^4$ 2x2 Matrices

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 ...
  • 8,521
7 votes
Accepted

roots of polynomials with small coefficients

Note that, if $D$ is invertible, the eigenvalues of $A$ and $DAD^{-1}$ are the same. You can avoid floating-point underflow when forming the matrix by scaling the companion matrix by a diagonal ...
  • 372
7 votes

Speeding up a linear transform using Python

Assuming that your kernel is somewhat smooth, use low-rank approximation. Here's a naive example: ...
  • 596
7 votes

Calculating partial trace of array in NumPy

I didn't follow exactly your notation so I can't say for sure what this would look like for your example, but two suggestions: if you really want your life made simple, check out ...
  • 191
7 votes

Numpy FFT gives me a pulse shorter than it should be. Not sure what I am doing wrong

Running your code, it seems like your pulse looks kinda like this: (sorry for not adding units to the plots, I used the same as you, i.e. t is in fs and w in rad/fs) So, the FWHM is not correct (...
7 votes

Is there an efficient way to form this block matrix with numpy or scipy?

The code proposed by the OP can indeed made be more efficient, mainly by noting the fact that to form the sequence $A^i B$, with $i=0\,\dots,N$ you do not have to compute $A^i$ at each step, but you ...
  • 3,809
7 votes
Accepted

Composite matrices in Numpy

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

Going to try to move some of my scipy/numpy calculation to a new GPU, how to avoid disappointing results?

I buy the wrong CUDA GPU and the speed up is minimal or nonexistent. It is highly unlikely that your choice of GPU will have a significant impact on your speed-up unless your model is very big. To a ...
  • 3,566
6 votes
Accepted

Python Vectorizing a Function Returning an Array

The problem is that np.cos(t) and np.sqrt(t) generate arrays with the length of t, whereas ...
  • 176
6 votes

What does Python offer for distributed/parallel/GPU computing?

Here's a few options that are relatively easy to work with: One node - multiprocessing is the most straightforward thing to do. multiprocessing.map works well for ...
6 votes

Complexity of matrix inversion in numpy

You should probably note that, buried deep inside the numpy source code (see https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src) the inv routine attempts to call the dgetrf ...
  • 2,199
6 votes
Accepted

Python: vectorizing a structured linear system solve

In your explanation, you solve the large problem using forward substitution. This implies that you are solving your large problem successively: you first need $x_{i-1}$ before you can solve for $x_{i}$...
  • 6,099
6 votes

Diagonalize a unitary matrix with orthogonal matrices using numpy

Have you tried the QZ decomposition on real(U) and imag(U)? In general it returns AA and <...
6 votes

Numerical derivative in python

The Savitzky-Golay filter uses a constant delta (the spacing of the samples,) and the default value of the delta in the filter implementation is 1, according to https://docs.scipy.org/doc/scipy-0.16.1/...
  • 71
6 votes

solve_ivp from scipy does not integrate the whole range of tspan

You can examine the sol object to see why the integration failed. It provides the message 'Required step size is less than spacing between numbers.' This usually ...
6 votes
Accepted

Calculate determinant of unitary matrices based on SVD implementation

If you are prepared to go digging around in the fortran code: The SVD algorithm consists of a few parts: Bidiagonalization (usually using Householder reflectors) Use QR shifts to reduce the ...
  • 1,254
5 votes

Python, numpy and complex functions (PDE's)

Caveat: I did not read beyond the statement So, what I did was to define functions for the potential and the second derivative, and use Euler's method. So there may be other issues with your code, ...
5 votes

Calculating partial trace of array in NumPy

I don't have enough reputation to comment, but I wanted to give a basic complement of gIS's answer. Simply put, if we have the decomposition over $V=V_1\otimes V_2$, which we represent with, say, the ...
  • 51
5 votes

Numerical stability in the product of many matrices

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 ...
  • 1,097
5 votes

How to Invert a Poorly Conditioned Matrix

There is no simple fix. For an ill-conditioned matrix $A$, the harm (loss of precision) is already done the moment you wrote those numbers in a numpy array, because that tiny $10^{-16}$ perturbation ...
5 votes
Accepted

Lanczos algorithm for finding top eigenvalues of a matrix sum

You could define a linear opearator and pass it to the function eigsh. Ideally, your matrices $L$ and $X$ are sparse so you can take advantage of the matrix-vector ...
  • 8,219
5 votes

Calculate determinant of unitary matrices based on SVD implementation

For part 1: To my knowledge, the answer is no. For part 2: This question had several good answers, all of which were negative (there isn't really a faster way). I don't believe there is meaningful new ...
5 votes
Accepted

Stochastic SIR using SDEint python package

This model is implemented using Julia's DifferentialEquations.jl in this tutorial. Here's a version of that code: ...
4 votes
Accepted

Performing a random walk on a lattice that traps the particles

This can be solved much more efficiently if you recast the problem as finding the mean first passage time (MFPT) of a Markov chain for a given starting state. You can easily represent the random walk ...
4 votes

Fast Automatic Differentiation for numpy?

Jax has the features you're looking for. See https://jax.readthedocs.io/en/latest/notebooks/quickstart.html
  • 149
4 votes

Moore-Penrose pseudoinverse of singular rank degenerate matrix

75k x 75k double-precision entries is 45 gigabytes. That fits in memory, but barely; you need to be careful. The linear algebra routines in most languages rely on Lapack as a backend, which is a ...

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