Questions tagged [numpy]
NumPy is the fundamental package for scientific computing with Python.
131
questions
30
votes
5answers
22k views
Permute a matrix in-place in numpy
I want to modify a dense square transition matrix in-place by changing the order of several of its rows and columns, using python's numpy library. Mathematically this corresponds to pre-multiplying ...
21
votes
1answer
2k views
How does the performance of Python/Numpy array operations scale with increasing array dimensions?
How do Python/Numpy arrays scale with increasing array dimensions?
This is based on some behaviour I noticed while benchmarking Python code for this question: How to express this complicated ...
15
votes
1answer
3k views
Why does SciPy eigsh() produce erroneous eigenvalues in case of harmonic oscillator?
I'm developing some larger code to perform eigenvalue computations of huge sparse matrices, in the context of computational physics. I test my routines against the simple harmonic oscillator in one ...
14
votes
3answers
2k views
How to express this complicated expression using numpy slices
I wish to implement the following expression in Python:
$$
x_i = \sum_{j=1}^{i-1}k_{i-j,j}a_{i-j}a_j,
$$
where $x$ and $y$ are numpy arrays of size $n$, and $k$ is a numpy array of size $n\times n$. ...
13
votes
2answers
9k views
Complexity of matrix inversion in numpy
I am solving differential equations that require to invert dense square matrices. This matrix inversion consumes the most of my computation time, so I was wondering if I am using the fastest algorithm ...
13
votes
1answer
1k views
Can an approximated Jacobian with finite differences cause instability in the Newton method?
I have implemented a backward-Euler solver in python 3 (using numpy). For my own convenience and as an exercise, I also wrote a small function that computes a finite difference approximation of the ...
10
votes
4answers
4k views
Memory efficient implementations of partial Singular Value Decompositions (SVD)
For model reduction, I want to compute the left singular vectors associated to the - say 20 - largest singular values of a matrix $A \in \mathbb R^{N,k}$, where $N\approx 10^6$ and $k\approx 10^3$. ...
9
votes
1answer
1k views
Numerical integration for modelling curve for superconductors (Python)
I am a physicist who is trying to model the current-voltage characteristics of a superconductor-superconductor junction.
The equation for this model is:
\begin{align}
I(V) = \frac{1}{eR_{\mathrm{n-n}...
8
votes
3answers
4k views
What does Python offer for distributed/parallel/GPU computing?
Using the SciPy/NumPy libraries, Python is a pretty cool and performing platform for scientific computing. I just wonder: When I have to go parallel (multi-thread, multi-core, multi-node, gpu), what ...
8
votes
2answers
631 views
Continuity of eigenvectors of parametric matrix
I have $n$-dimensional matrices $\mathrm{\hat{H}}(\vec{k})$ depending on vector parameter $\vec{k}$.
Now, eigenvalue routines return eigenvalues in no particular order (they are usually sorted), but ...
8
votes
0answers
590 views
Fast Automatic Differentiation for numpy?
I would like to use automatic differentiation to calculate gradients to function written in numpy.
I've come across a number of packages, including
autograd
tangent
chainer
But none of them seem ...
7
votes
3answers
12k views
plotting discontinuous functions
I need help plotting the Heaviside function:
Real analysis often involves constructing bizarre functions which are intuitively correct, but ultimately wrong. See the great book Counterexamples in ...
7
votes
0answers
166 views
fastest way to compute many small dot products
I have two n-by-3 blocks contiguous in memory ("n vectors of length 3") and I'd like to compute the dot product between each of the rows as fast as possible. In numpy, using ...
6
votes
3answers
4k views
Calculating partial trace of array in NumPy
A simulation I'm doing requires me to calculate the partial trace of a large density matrix. I am trying to calculate it using tools from numpy, but my code seems to be having some problems. For ...
6
votes
2answers
2k views
Is there an efficient way to form this block matrix with numpy or scipy?
Is there an efficient way to form this block matrix with numpy or scipy?
$$
\left[
\begin{array}{cccc}
\mathbf{B} & \mathbf{0} & \cdots & \mathbf{0}\\
\mathbf{AB} & \mathbf{B} & \...
6
votes
1answer
385 views
roots of polynomials with small coefficients
I would like to compute the roots of a polynomial with exponentially small coefficients.
$$ \sum_{n=0}^N a_n \frac{z^n}{\sqrt{n!}} \tag{$\ast$}$$
where $a_n$ are Normal random variables with mean $0$...
6
votes
2answers
4k views
Using numpy arrays in Paraview programmable filter
How can I access a field in Paraview's programmable filter as a numpy array?
I want to:
Import an existing field as a numpy array
Create a similar array for output
Register it as a new field for ...
6
votes
1answer
400 views
Calculate large and small frequency separation for the Sun
I want to determine the big and small frequency seperation from timeseries data for the sun. An excerpt of the data (timeseries and power series) is plotted below.
The power series is calculated in ...
5
votes
1answer
2k views
Alternatives to numpy.einsum
Given an $n_1 \times \cdots \times n_k \times g \times g$ tensor $A$ (i.e. a collection of $g \times g$ matrices) and an $n_1 \times \cdots \times n_k \times g$ tensor $b$ (i.e. a collection of ...
5
votes
2answers
387 views
Speeding up a linear transform using Python
In an optical wavefront propagation problem, I need to do excessive Fourier-type computations:
...
5
votes
2answers
404 views
Python: vectorizing a structured linear system solve
Overview
I am looking for a way to solve a structured linear system in Python without using a for loop (preferably using vectorization, if possible).
Background
Consider the following linear system:...
5
votes
1answer
2k views
Sign differences in spectral decomposition in NumPy
I am trying to understand an example from a book, but I seem to get different answers depending on which spectral decomposition function I use in NumPy.
I am trying to find a spectral decomposition $...
5
votes
1answer
1k views
Finite Volume Implementation
I am trying to implement a simple finite volume method solver. I had a class on FVM a while back, but am still aware of the principal concepts. But implementing the FVM for non-cartesian or 1D meshes ...
5
votes
1answer
3k views
Logistic regression with Python
I am trying to code up logistic regression in Python using the SciPy fmin_bfgs function, but am running into some issues. I wrote functions for the logistic (...
5
votes
1answer
92 views
How do I globally change the precision of a piece of code in Python to debug it?
I am solving a system of non-linear equations using the Newton-Raphson method in Python. This involves using the solve(Ax,b) function (...
5
votes
0answers
1k views
Python - calculation time derivative and laplacien by finite differences
I would like to determine a temporal derivative and a Laplacian by the finite differences method to solve the heat equation in a 1-dimensional case. My aim is to get the sources term that is why I ...
4
votes
3answers
6k views
Can numpy.linalg.solve use back substitution when possible?
The question is if Python Numpy library can use back subsitution to solve Ax=b if possible, that is, if A is lower triangular? Do numerical linear algebra packages do this? I would think Numpy would ...
4
votes
2answers
217 views
Fastest Way to Mutiply $10^4$ 2x2 Matrices
In a code that I work with (written in python, but also tagging as matlab because numpy is so close and I could use it if need be), we use a transfer matrix method to compute the properties of a ...
4
votes
1answer
653 views
Poor SVD reconstruction of singular matrix
I am trying to calculate the singular value decomposition of this matrix using numpy.linalg.svd .
However, reconstructing the matrix from the SVD gives a poor ...
4
votes
4answers
10k views
LCM builtin in Python / Numpy
I can write a function to find LCM (lowest common multiple) of an array of integers, but I thought it must have been implemented in numpy or scipy and was expecting something like ...
4
votes
2answers
9k views
Consistent handling of division by zero in numpy array
I want to populate a numpy array with values from the smooth bump function
f(x) = exp ( - 1 / (1 - x^2) ) if |x| < 1, f(x) = 0 otherwise
Currently I ...
4
votes
1answer
393 views
Diagonalize a unitary matrix with orthogonal matrices using numpy
An important component of the Cartan KAK decomposition for 2 qubit operations is to diagonalize a 4x4 unitary matrix using orthogonal (not unitary, purely real orthogonal) matrices. That is to say, ...
4
votes
2answers
355 views
Moore-Penrose pseudoinverse of singular rank degenerate matrix
I am trying to attain the Moore-Penrose pseudoinverse of a very large, very sparse, rank-degenerate, singular, and square matrix. ($75000 \times 75000$, near rank). The matrix is a graph Laplacian and ...
4
votes
1answer
1k views
Python syntax for MATLAB/Octave colon operator a:dx:b
I am trying to rewrite some MATLAB/Octave code in Python, and I don't know what would be the nicest or most intuitive way of writing
...
4
votes
1answer
64 views
Backwards Difference Implicit Method for Nonlinear Parabolic PDE in Python
Original Stack Overflow Question: https://stackoverflow.com/questions/65683788/indexerror-index-31-is-out-of-bounds-for-axis-1-with-size-31?noredirect=1#comment116218335_65683788
PDE: u_t = u_xx + u(...
4
votes
1answer
655 views
Chebyshev approximation by projection vs interpolation
Suppose we want to approximate a function $f: [a, b] \rightarrow \Re$ with a Chebyshev series:
$$ f(x) \approx \sum_{k=0}^n c_k \, T_k\left( \frac{2x-b-a}{b-a} \right) $$
where $T_k(x) = \cos(k\, \...
4
votes
2answers
106 views
How to read the number of periods of this complicated graph?
I have two data sets that are quasi-periodic. They have the same period and can be seen clearly by eye. For example when $x\in(100,200)$, both of them have about 32 periods.
Below is a zoom-in of the ...
4
votes
1answer
2k views
Analytic solution 2D scalar wave equation in cylindrical coordinates numerical implementation
I am trying to compare my finite difference's solution of the scalar (or simple acoustic) wave equation with an analytic solution.
For that purpose I am using the following analytic solution ...
3
votes
1answer
230 views
High frequency noise at solving diffusion equation
I'm trying to simulate a simple diffusion based on Fick's 2nd law.
...
3
votes
3answers
5k views
Runge-Kutta Simulation For Projectile Motion With Drag
I am attempting to simulate projectile flight with drag. However, with a timestep of 0.1 seconds, I am consistently getting an error of ~0.1-1%.
...
3
votes
3answers
2k views
Caculating the mean of vector accurately
I am having trouble with calculating a mean of vector with sufficient accuracy. My current solution which works but it quite slow and has unpredictable performance:
mean_sum = mean = math.fsum(values)...
3
votes
1answer
991 views
Numerical derivative in python
I am trying to take the numerical derivative of a dataset. My first attempt was to use the gradient function from numpy but in that case the graph of the derivative ...
3
votes
1answer
6k views
Python Vectorizing a Function Returning an Array
I have the following function that has been vectorized so that for every element in input array t, an array is output:
...
3
votes
1answer
118 views
Fast nonzero indices per row/column for (sparse) 2D numpy array
I am looking for the fastest way to obtain a list of the nonzero indices of a 2D array per row and per column. The following is a working piece of code:
...
3
votes
1answer
59 views
why am I not getting a staircase for the rotation number?
I'm trying to understand the staircase map. Look at this map from the circle to itself:
$$ x \stackrel{F}{\mapsto} \big[\omega + x + \tfrac{\epsilon}{2\pi} \sin (2\pi x) \big] \pmod 1 $$
Such a map ...
3
votes
1answer
665 views
My calculated laser pulse duration is too large. Where am I wrong?
I am currently writing a small Python script to estimate the pulse duration from the optical spectrum.
At the end, the idea is to observe the effects of the spectral phase on the pulse duration and ...
3
votes
1answer
940 views
32bit/64bit issue when working with Numpy and petsc4py
When indexing PETSc.Mat A with an array c ( numpy.ndarray with ...
3
votes
1answer
1k views
integral including a vector
I need to calculate the integral of this function
def f(z):
return ((1-2*z)*np.exp(-d/z))/(((1-z)**(2+d))*(z**(2-d)))
Here d is a constant. I am using this ...
3
votes
1answer
156 views
Riemann surfaces: computing $f(z) = \int_0^z \frac{dx}{\sqrt{P(x)}}$
I am trying to validate that the Schwartz-Christoffel mapping does indeed take the upper-half plane $\mathbb{H} = \{ z: \mathrm{Re}(z) > 0\}$ to a polygon.
This involves integrals of functions ...
3
votes
0answers
36 views
Difference between wave vector and density matrix in numerical calculation of Schrödinger equation
I solved Schrödinger equation for a following tow-level time-dependent Hamiltonian numerically in two ways:
...
