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 ...
- 12k
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^* - ...
- 9,666
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 (...
- 256
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.
- 9,666
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 <...
- 9,666
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 ...
- 989
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, ...
- 16.3k
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 ...
- 9,666
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 ...
- 1,357
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:
...
- 12k
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 ...
- 3,954
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 ...
- 9,666
Only top scored, non community-wiki answers of a minimum length are eligible