Questions tagged [numpy]
NumPy is the fundamental package for scientific computing with Python.
159
questions
31
votes
5
answers
28k
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
1
answer
3k
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 ...
17
votes
2
answers
13k
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 ...
16
votes
1
answer
4k
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
3
answers
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
1
answer
1k
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 ...
13
votes
1
answer
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 ...
10
votes
4
answers
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
2
answers
822
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 ...
9
votes
1
answer
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
3
answers
10k
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 ...
8
votes
3
answers
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
3
answers
18k
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
0
answers
372
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
3
answers
9k
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 ...
6
votes
1
answer
3k
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
2
answers
5k
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
1
answer
474
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
2
answers
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
1
answer
448
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
1
answer
2k
views
Rational functions in Python
I would like use Python to verify the following identities:
$ \frac{1}{1-z} = 1 + z + z^2 + z^3 + \dots $
$ \frac{1}{1-z - z^2} = 1 + z + 2z^2 + 3z^3 + \dots $
$ q \prod_{n \geq 1} (1 - q^n)^{24} = ...
5
votes
2
answers
855
views
Speeding up a linear transform using Python
In an optical wavefront propagation problem, I need to do excessive Fourier-type computations:
...
5
votes
2
answers
519
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
1
answer
589
views
Stochastic SIR using SDEint python package
I want to use the SDEint package to give a numerical solution (plot) of the following stochastic SIR model. Namely, a system of SDEs.
$$\begin{cases}
dS = -\beta SIdt - \sigma SIdW \\
dI = (\beta SI -...
5
votes
1
answer
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
1
answer
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
1
answer
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
1
answer
152
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
0
answers
98
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 ...
5
votes
0
answers
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
2
answers
620
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
1
answer
895
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
4
answers
11k
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
2
answers
12k
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
1
answer
357
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....
4
votes
1
answer
977
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
2
answers
533
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
1
answer
2k
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
1
answer
665
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
1
answer
920
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
2
answers
109
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
0
answers
159
views
How amenable is this 2D Frenkel–Kontorova-like energy minimization problem in Python to the use of a modest PC + GPU? (Heavy reliance on indexing)
@Richard's answer to Going to try to move some of my scipy/numpy calculation to a new GPU, how to avoid disappointing results? is quite helpful, and as promised I've added a simple running example ...
4
votes
0
answers
137
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
1
answer
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
1
answer
169
views
Time and memory required to diagonalize a 18000 by 18000 matrix using numpy in python
Can someone give an estimate of the Time and memory required to diagonalize a 20000 by 20000 complex hermitian matrix using numpy in python ?
3
votes
1
answer
239
views
High frequency noise at solving diffusion equation
I'm trying to simulate a simple diffusion based on Fick's 2nd law.
...
3
votes
3
answers
8k
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
3
answers
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
1
answer
3k
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
1
answer
8k
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:
...