22
(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 claimed to be faster, just that they have better asymptotic complexity (meaning the required number of elementary operations grows more slowly). In fact, in ...
17
Here is the Numba solution. On my machine the Numba version is >1000x faster than the python version without the decorator (for a 200x200 matrix, 'k' and 200-length vector 'a'). You can also use the @autojit decorator which adds about 10 microseconds per call so that the same code will work with multiple types.
from numba import jit, autojit
@jit('f8[:]...
12
The degeneracy of some eigenvalues looks to me like the hallmark of the breakdown of the Lanczos algorithm. The Lanczos algorithm is one of the more commonly used methods to approximate the eigenvalues and eigenvectors of Hermitian matrices; it's what scipy.eigsh() uses, through a call to the ARPACK library.
In exact arithmetic, the Lanczos algorithm ...
9
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 that NumPy will adopt for this operator is "stacked matrix multiplication," i.e. exactly what you want.
There may be bits and pieces missing, but some work is already in, see the relevant PR, and ...
9
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^* - A \|}{\|A\|} = O(u)$, where the norms are Euclidean norms, $u$ is the machine precision, and "$O(u)$" means that the first-order term in $u$ is bounded by a ...
8
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 ensure that your timestep is small enough compared to the grid spacing.
In the case of the heat equation (Fick's 2nd law is a heat equation), your CFL number $C$ ...
8
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 (as mentioned by @Mauro Vanzetto) as well as particular libraries/packages.
I, personally, like having np.sqrt or std::cout (as opposed to just sqrt and cout) ...
7
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 matrix $D$ such that $D_{ii} = 1/\sqrt{(n-i)!}$. For the polynomial $x^4 + a x^3/\sqrt{1!} + b x^2/\sqrt{2!} + c x/\sqrt{3!} + d/\sqrt{4!}$, this gives a modified ...
7
Assuming that your kernel is somewhat smooth, use low-rank approximation.
Here's a naive example:
import numpy as np
N=2000
input=np.random.random(N)
x=np.linspace(-1,1,N)
y=np.linspace(-2,2,N)
X,Y=np.meshgrid(x,y,sparse=True)
A = np.exp(1j*2*np.pi*X*Y)
output = np.dot(A, input)
U,S,V = np.linalg.svd(A)
# find truncation rank for given tolerance
k = ...
7
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 (should be 4 fs) and the angular carrier/central frequency is messed up (should be ca. 2.3 / ca. 0.37 per fs). I think there is two things not quite right here:
1) ...
7
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 can exploit the fact that $A^i B = A\,(A^{i-1}B)$, reusing the result of the previous step.
My proposed implementation is
import numpy as np
N = 3
M = 1
A = ...
7
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 numpy:
consider np.array vs np.matrix, it might happen that np.matrix is faster than np.array matrix-matrix product (it is unclear what you are using now, and ...
7
They are commonly called block matrices. You can create them with hstack, vstack, and block.
6
How accurately? Is there a reason you're not just using numpy.mean? Is that not sufficient? If you need to compensate for floating-point error, you can try using Kahan summation. Here's some pseudocode (PDF).
For a review of floating-point summation techniques see: N.J. Higham, SIAM, 1993.
6
No. The numpy.linalg.solve method uses LAPACK's DGESV, which is a general linear equation solver driver. If you know that your matrix is triangular, you should use a driver specialized for that matrix structure.
scipy.linalg.solve does something similar.
MATLAB detects triangularity in a solve if you use the backslash operator; see this page for pseudocode....
6
If you only want a few singular values/vectors, ARPACK should do the trick. The SVD docs aren't great, and this distribution is more up to date.
EDIT: If you want to do this in python, SciPy has a wrapper. Since your matrix is dense, you could try the block sparse row (BSR) format.
6
Warning: The below example works properly, but using the full set of parameters suggested at the post end exposes a bug, or at least an "undocumented feature" in the numpy.take() function. See comments below for details. Bug report filed.
You can do this in-place with numpy's take() function, but it requires a bit of hoop jumping.
Here is an example of ...
6
def f(r):
return 2 * np.clip(r, 0, 1) / (r + 1)
6
Looking at the information of nympy.linalg.solve for dense matrices, it seems that they are calling LAPACK subroutine gesv, which perform the LU factorization of your matrix (without checking if the matrix is already lower triangular) and then solves the system. So the answer is NO.
Otherwise, it makes sense. If you do not have an easy (cheap) way to verify ...
6
The problem is that np.cos(t) and np.sqrt(t) generate arrays with the length of t, whereas the second row ([0,1]) maintains the same size.
To use np.vectorize with your function, you have to define the output type, and np.vectorize
isn't really meant as a decorator except for the simplest cases.
In this way however you can generate the function with the ...
6
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 an embarrasingly parallel problem.
GPU - pyCUDA allows you to do some GPU level programming with just python. I haven't toyed with it that much.
MPI - mpi4py- Exposes MPI interface at the python ...
6
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}$. This means that you need to loop, there is no way to avoid this.
However, for each subproblem, you use the spsolve routine, dedicated to solving sparse ...
6
Have you tried the QZ decomposition on real(U) and imag(U)? In general it returns AA and BB upper triangular rather than diagonal, but I wonder if a stroke of luck happens here (exactly like the Schur decomposition can be proved to always return a diagonal $T$ when run on a normal matrix, rather than an upper triangular one).
(EDIT: yes, on second thought ...
6
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/reference/generated/scipy.signal.savgol_filter.html. Your data set has irregular deltas, not 1, so the result from the Savgol filter is incorrect.
5
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 values for L1 and L2 cache sizes of a (a bit old) AMD CPU. To be certain, I did my own benchmark:
def timeit(size):
t0 = time.time()
for _ in xrange(10):
...
5
The "fast" version is (damped/over-relaxed) Jacobi instead of SOR, which is a multiplicative method. Even with damping factor of 1.0, Jacobi is not guaranteed to converge, as you would see if you applied it to a 1D problem, or added extreme anisotropy to your multi-dimensional problem.
Your SOR ("slow") implementation is slow because it is written in pure ...
5
Why do you want to generate DivisionByZero exceptions?
I would use masked arrays:
import numpy as np
x= np.linspace(-1.1,1.1,300)
masked_idx = (np.abs(x)>1)
masked_x = np.ma.array(x,mask=idx)
def f(x):
return np.exp(-1.0/(1.0-x**2))
masked_f = f(masked_x)
plot(masked_x,masked_f) # in IPython / pyplot
If you want, you can do the masking in ...
5
I think you are looking for capabilities not available in Python itself nor Numpy, but that are available in SymPy. Sympy is a computer algebra package developed in Python. You can try it "live" in its on-line shell. Documentation on series development can be found here.
5
Let $w_k$ be the k-th column of W (the kth eigen-vector) and $v_k$ be the k-th element of v (the kth eigen-value). Then we can write:
$$ G = \sum_k w_k v_k w_k^T $$
which, element-wise, is equivalent to:
$$ G_{i,j} = \sum_k v_k w_{i,k} w_{j,k} $$
where the indices are column-major.
Since we always get $w_k$ twice, an overall sign on a column of $W$ cancels.
...
5
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, but this one is already fundamental.
You are trying to solve a wave equation with Euler's method. This is numerically unstable. To see why, note that the ...
Only top voted, non community-wiki answers of a minimum length are eligible
