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 ...
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
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 ...
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
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....
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
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 qutip, the Quantum Toolbox in Python. It has classes specifically for quantum operators and state vectors, supports tensor product of operators in different spaces and keeps track ...
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
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
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
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 function from your system LAPACK package, which then performs an LU decomposition of your original matrix. This is morally equivalent to Gaussian elimination, but ...
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.
6
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 bidiagonal to diagonal (usually using Givens rotations)
Change the sign of any negative diagonal elements and reorder the singular values so the are decreasing.
...
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
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
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
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
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 ...
5
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 after some number of multiplications. You can do this by just taking the QR decomposition of your matrix after some number of products and taking the orthogonal ...
5
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 from the exact non-representable values is already harmful.
You could increase your working precision; but at that point the question is if your matrix entries $...
Only top voted, non community-wiki answers of a minimum length are eligible