Questions tagged [numpy]

NumPy is the fundamental package for scientific computing with Python.

Filter by
Sorted by
Tagged with
30
votes
5answers
26k 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 ...
16
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
11k 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
2k 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 ...
11
votes
1answer
708 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 ...
10
votes
4answers
5k 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
709 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 ...
7
votes
3answers
6k 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 ...
7
votes
3answers
14k 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
182 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
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 ...
6
votes
2answers
4k 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
425 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
427 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
2answers
569 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
443 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
107 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
7k 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
445 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
779 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
11k 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
680 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
425 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
276 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
756 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
107 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
0answers
40 views

How to troubleshoot numerical instability using finite difference for steady-state non-linear heat conduction equation

I have a problem which I believe is numerical instability when trying to solve a heat conduction equation using finite difference. The short version is that when the parameter $I=80.3$ I get the blue ...
4
votes
0answers
57 views

Solving multiple linear regression in parallel

I am working on a problem where I need to solve approximately 500 Million Linear Regressions (OLS). What would be the most efficient way to do this (e.g. using GPU or a some framework that can do this ...
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
237 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
6k 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
3k 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
2k 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
7k 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
109 views

Lanczos algorithm for finding top eigenvalues of a matrix sum

I am trying to find the top k leading eigenvalues of a NumPy matrix (using python dot product notation) L@L + a*X@X.T, where $L$ ...
3
votes
1answer
328 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
131 views

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

update: I've refactored the question based on helpful advice in the linked meta. I'm a heavy user of Python's NumPy and SciPy (and not much else) and for years I could run anything I need on my laptop....
3
votes
1answer
40 views

Forming a particular (averaged) block matrix with numpy

Say I have a set of $n \times n$ matrices $A_1, ..., A_m$ as numpy arrays. I'd like to create the block matrix defined below. I'm looking for a clean, elegant, and easy-to-interpret way of doing this ...
3
votes
1answer
63 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 ...