60
Except for code which does a significant number of floating-point operations on data that are held in cache, most floating-point intensive code is performance limited by memory bandwidth and cache capacity rather than by flops.
$v$ and the products $Av$ and $Bv$ are all vectors of length 2000 (16K bytes in double precision), which will easily fit into a ...
27
Edit: This is now in SymPy
$ isympy
In [1]: A = MatrixSymbol('A', n, n)
In [2]: B = MatrixSymbol('B', n, n)
In [3]: context = Q.symmetric(A) & Q.positive_definite(A) & Q.orthogonal(B)
In [4]: ask(Q.symmetric(B*A*B.T) & Q.positive_definite(B*A*B.T), context)
Out[4]: True
Older answer that shows other work
So after looking into this for a while ...
20
Your matrix is of size 15,000 x 15,000, so you have 225M elements in the matrix. This makes
for roughly 2GB of memory. This is much more than the cache size of your processor, so it has to be loaded completely from main memory in every matrix multiplication, making for approximately 100GB of data transfers, plus what you need for the source and destination ...
19
See https://math.stackexchange.com/questions/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (sorry, I would have put that in a comment but I've registered just to post this so I can't post comments yet).
But since I'm writing it as an answer, I'll also write the method:
$$E=\frac{m_{00}+m_{11}}{2}; F=\frac{m_{00}-m_{11}}{2}; G=\...
19
Your code is limited by memory bandwidth. For trivial math, it's often better to count memory accesses rather than flops. You'll get the following table:
operation memory reads/writes
matrix + matrix 3n²
matrix * vector 2n²+n (if vector is not cached)
matrix * vector n²+2n (if vector is only read once)
vector + vector ...
18
The column major layout is the scheme used by Fortran and that's why it's used in LAPACK and other libraries.
In general it is much more efficient in terms of memory bandwidth usage and cache performance to access the elements of an array in the order in which they're laid out in memory. Depending on how your matrices are stored, you'll want to pick ...
17
Matlab interprets sequences of multiplications and/or divisions from left to right. Hence $A*B*C*v$ is much more expensive than $A*(B*(C*v))$, as you have two matrix products and one matrix-vecor product in place of three matrix-vector products.
On the other hand, $A*(B*(C*v))$ should be slightly faster than if you save the intermediates in separate ...
17
The first thing is to recognize that you can do this using BLAS. If you data matrix is $X = [x_1 x_2 x_3 ...] \in \mathbb{R}^{m\times n}$ (each $x$ is a column vector corresponding to one measurement; rows are trials), then you can write the covariance as:
$$ C_{ij} = E[x_i,x_j] - E[x_i] E[x_j] = \frac{1}{n} \sum_k x_{ik} x_{jk} - \frac{1}{n^2} \left(\sum_k ...
14
For starters, I wouldn't use intermediate variables, but brackets. Unless, of course, you're interested in the intermediate results, but I'm guessing not.
I tried the following in Matlab:
>> N = 500;
>> A = rand(N); B = rand(N); C = rand(N); v = rand(N,1);
>> tic, for k=1:100, A*B*C*v; end; ...
14
Since the matrices are so small, all of the cost is going to be in call overhead. If you will do the transformation many times, it will be faster to precompute D=A*B*C once and then for each vector apply v_f=D*v_i. You could also consider bringing this out to a mex file.
14
You can just simulate the matrix-matrix product by forming the product of the two sparsity patterns -- i.e., you consider the sparsity pattern (that is stored in separate arrays in CSR format) as a matrix that contains either a zero or a one in each entry. Performing this simulated product only requires you to form the and operation on these zeros and ones ...
14
The Cholesky factorization $C=R^TR$ leads to a Cholesky-like factorization of the inverse $C^{-1}=SS^T$ with the upper triangular matrix $S=R^{-1}$.
In practice, is best to keep the inverse factored. If $R$ is sparse then it is usually even better to keep $S$ implicit, as matrix-vector products $y=C^{-1}x$ can be computed by solving the two triangular ...
14
Matrix Market is a terrible format for reading in parallel, therefore it is better to preprocess to a better parallel format. Your matrix size is extremely small so performance is not an issue, but the easiest and most general thing is to use Python or Matlab/Octave to write the Matrix Market file in PETSc binary format, which can be read efficiently in ...
14
The most obvious thing you can do is to precompute
[L,U] = lu(A) ~ O(n^3)
Then you just compute
x = U \ (L \ b) ~ O(2 n^2)
This would reduce the cost enormously and make it faster. Accuracy would be the same.
14
Here is R1, as computed in MATLAB:
1.0e+07 *
-7.382605957465515 -9.599867106092937 -2.830412177259742 -0.000000000002830 -0.000000000002830
-1.230434326244253 -1.599977851015490 -0.471735362876624 -0.000000000000472 -0.000000000000472
3.691302978732758 4.799933553046468 1.415206088629871 0.000000000001415 0.000000000001415
-5....
14
First, see Mark L. Stone's answers, which is completely correct. Second, realize that this is the reason why people told you to use relative errors in your numerical analysis class. :)
Third, the real question here is why the results do not coincide exactly, since both languages call some BLAS library functions for their computations. There are several very ...
13
I actually wrote the original code in Matlab for A*B, both A and B sparse. Pre-allocation of space for the result was indeed the interesting part. We observed what Godric points out -- that knowing the number of nonzeros in AB is as costly as computing AB.
We did the initial implementaion of sparse Matlab around 1990, before the Edith Cohen paper that ...
13
Apart from the (now classical) Golub-Reinsch paper Brian notes in his answer (I have linked to the Handbook version of the paper), as well as the (also now classical) predecessor paper of Golub-Kahan, there have been a number of important developments in computing the SVD since then. First, I have to summarize how the usual method works.
The idea in ...
13
In fundamental C++, I find the problem here is that C++ will allocate a new object of cx_mat to store evolutionMatrix*stateMatrix, and then copy the new object to stateMatrix with operator=().
I think you're right that it's creating temporaries, which is too slow, but I think the reason for why it's doing that is wrong.
Armadillo, like any good C++ linear ...
13
If the solution of $Ax=b$ is unstable, the matrix is very ill-conditioned (i.e., has a very large condition number), and (paraphrasing Lanczos) no amount of mathematical trickery can make it stable. The best you can hope for is to solve a different problem that is a) stable and b) gives you a solution that is sufficiently close; this is called regularization....
12
In gerschgorin's theorem, the diagonal entries $A_{ii}$ of the matrix are the eigenvalue estimates, and the radii $r_i$ of the Gerschgorin disks are corresponding error bounds. Thus $\min_i A_{ii}-r_i$ is a lower bound on the eigenvalues, and $\max_i A_{ii}+r_i$ is an upper bound.
Note that these bounds are generally poor unless the off-diagonal entries are ...
12
When computing the eigenvalues of the symmetric matrix $M\in\mathbb{R}^{n\times n}$ the best you can do with Householder reflector is drive $M$ to a tridiagonal form. As was mentioned in a previous answer because $M$ is symmetric there is an orthogonal similarity transformation which results in a diagonal matrix, i.e., $D=S^TMS$. It would be convenient if we ...
12
The LU factors of a sparse matrix are at least somewhat sparse. The $Q$ matrix in QR can also somewhat preserve sparsity, and is typically used when the matrix is very long and skinny. The SVD of a sparse matrix will almost always have fully dense $U$ and $V$ factors, so it destroys any reason to perform the computations treating the matrix sparsely.
12
The property follows from the property of the corresponding (weak form of the) partial differential equation; this is one of the advantages of finite element methods compared to, e.g., finite difference methods.
To see that, first recall that the finite element method starts from the weak form of the Poisson equation (I'm assuming Dirichlet boundary ...
12
Here is a quick example which is very practical related to memory usage in PDEs. When one discretizes the Laplace operator, $\Delta u$, for example, in the Heat Equation
$$ u_t = \Delta u + f(t,u) .$$
To solve it numerically, one ends up with sparse matrices $A$, and a method of lines discretization then solve
$$ u_t = Au + f(t,u) $$
The canonical 1D ...
11
My previous answer recommended Dixon's 1983 paper, "Estimating extremal eigenvalues and condition numbers of matrices". It essentially boils down to a modest number of matrix-vector multiplications and solves against Gaussian random vectors and is essentially the power algorithm coupled with a priori error bounds which are not dependent on the spectrum of ...
11
As the Comments to other Answers clarify, the real issue here is not a shortcoming of Householder matrices but rather a question as to why iterative rather than direct ("closed-form") methods are used to diagonalize (real) symmetric matrices (via orthogonal similarity).
Indeed any orthogonal matrix can be expressed as a product of Householder matrices, so ...
11
We have the matrix Laplacian matrix $G=A^TA$ which has a set of eigenvalues $\lambda_0\leq\lambda_1\leq\ldots\leq \lambda_n$ for $G\in\mathbb{R}^{n\times n}$ where we always know $\lambda_0 = 0$. Thus the Laplacian matrix is always symmetric positive semi-definite. Because the matrix $G$ is not symmetric positive definite we have to be careful when we ...
11
In exact arithmetic you shouldn't need to reorthogonalize regularly, but practically you do. Your u1 and u2 are close to (but not exactly) the true eigenvectors, so your initial deflation almost (but not entirely) removed the true eigenvectors from u3. The tiny components you left behind will be amplified by repeated multiplication by A, you will need to ...
Only top voted, non community-wiki answers of a minimum length are eligible
