Hot answers tagged

25 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 ...
Christian Clason's user avatar
15 votes

Time and memory required to diagonalize a 18000 by 18000 matrix using numpy in python

A 20000 by 20000 double-precision complex matrix requires $20000 \times 20000 \times 8 \times 2=6.4 \mbox{gigabytes}$ of RAM. The LAPACK routines ZHEEV that will do the work for you will store the ...
Brian Borchers's user avatar
11 votes

What algorithm(s) do numpy and scipy use to calculate matrix inverses?

Documentation to numpy.linalg.inv and scipy.linalg.inv does not mention the algorithm used. Judging from the source, ...
Vladimir Lysikov's user avatar
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^* - ...
Federico Poloni's 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 (...
Anton Menshov's user avatar
  • 8,652
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 ...
Anton Menshov's user avatar
  • 8,652
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 ...
Okarin's user avatar
  • 191
7 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 ...
origimbo's user avatar
  • 2,239
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: ...
cfh's user avatar
  • 586
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 (...
matthiaw91's 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 ...
Stefano M's user avatar
  • 3,839
7 votes
Accepted

Composite matrices in Numpy

They are commonly called block matrices. You can create them with hstack, vstack, and block.
Federico Poloni's user avatar
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 ...
Richard's user avatar
  • 3,961
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 ...
agold's user avatar
  • 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 ...
Danielle's user avatar
6 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 ...
Cris's user avatar
  • 61
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}$...
GertVdE's user avatar
  • 6,179
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 <...
Federico Poloni's 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/...
KJ Nam's user avatar
  • 81
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 ...
Steven Roberts's user avatar
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 ...
Thijs Steel's user avatar
  • 1,558
6 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: ...
Chris Rackauckas's user avatar
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 ...
whpowell96's user avatar
  • 2,259
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 ...
Federico Poloni's 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 ...
nicoguaro's user avatar
  • 8,490
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 ...
Juan M. Bello-Rivas's user avatar
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
Brannon's user avatar
  • 148
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 ...
Federico Poloni's user avatar
4 votes
Accepted

Implementing structured grid boundary conditions using NumPy arrays?

The Numpy function $roll$ performs periodic shift of an array. Using it, the explicit time step for your PDE in a periodic domain can be simply implemented like this: ...
Maxim Umansky's user avatar

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