Hot answers tagged

23 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 ...
user avatar
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 ...
user avatar
  • 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^* - ...
user avatar
8 votes
Accepted

High frequency noise at solving diffusion equation

Solving partial differential equations with explicit timestepping methods relies on meeting a certain CFL condition for stability. Since you are using a forward Euler timestepping scheme, you must ...
user avatar
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 (...
user avatar
  • 8,382
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 ...
user avatar
  • 8,382
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 ...
user avatar
  • 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: ...
user avatar
  • 586
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 ...
user avatar
  • 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 (...
user avatar
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 ...
user avatar
  • 3,789
7 votes
Accepted

Composite matrices in Numpy

They are commonly called block matrices. You can create them with hstack, vstack, and block.
user avatar
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 ...
user avatar
  • 176
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 ...
user avatar
  • 2,199
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 ...
user avatar
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}$...
user avatar
  • 6,066
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 <...
user avatar
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/...
user avatar
  • 71
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 ...
user avatar
  • 1,189
6 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 ...
user avatar
  • 3,111
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, ...
user avatar
5 votes
Accepted

How does the performance of Python/Numpy array operations scale with increasing array dimensions?

I don't know how this benchmark was done, but probably with floating point numbers, that in Python default to doubles. The sizes correspond, respectively, to $4$ and $16 kB$. Those are reasonable ...
user avatar
  • 320
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 ...
user avatar
  • 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 ...
user avatar
  • 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 ...
user avatar
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 ...
user avatar
  • 8,006
5 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 ...
user avatar
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 ...
user avatar
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: ...
user avatar
4 votes

plotting discontinuous functions

As long as you know the exact positions of the discontinuities, you just have to set the jump positions to nan in x, ...
user avatar
  • 1,135

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