Hot answers tagged

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 ...


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 ...


9

Given $A \in {\bf S}^n$ (a positive definite matrix) with eigenvalues $\lambda_1 \leq \lambda_2 \leq \ldots \leq \lambda_n $, then: $\displaystyle f_k(A)=\sum_{i=1}^{k} \lambda_i$ is concave. Why? $$f_k(A) = \inf \left\{ {\bf tr}(V^T A V) | V \in {\bf R}^{n \times k}, V^T V = I \right\}$$ This follows from the Poincare separation theorem (see e.g. Horn ...


8

No. A reordering of a matrix $\boldsymbol{A}$ is equivalent to conjugation by a permutation matrix $\boldsymbol{P}$. In other words, the reordered matrix can be written as $\boldsymbol{A}_r = \boldsymbol{P}\boldsymbol{A}\boldsymbol{P}^*$ and $\boldsymbol{A}_r^{-1} = \boldsymbol{P}\boldsymbol{A}^{-1}\boldsymbol{P}^*$. Since the spectral norm is unitarily ...


8

Let $P$ be the anti-diagonal permutation matrix, $$P = \begin{bmatrix} & & & 1 \\ & & 1 \\ & 1 \\ 1 \end{bmatrix}$$ so that $PAP$ is the version of $A$ with rows and columns reversed. The first $P$ swaps the rows, and the last $P$ swaps the columns. We have the Cholesky decomposition $$PAP=LL^T$$ which implies $$A = (PLP)(PLP)^T,$$ ...


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

It's the blocked variant of Householder-QR that is driving this design. If you look in Golub and Van Loan's book (Ch 5.2 or so) they talk about how k-iterations of the algorithm can be blocked together by accumulating the individual reflectors into a rank-k reflector of the form $\mathbf I + \mathbf W \mathbf Y^{\mathrm T}$, where both $\mathbf W$ and $\...


6

As suggested in the comments, the eigendecomposition $\mathbf J = \mathbf Q \mathbf D \mathbf Q^{T}$ can be used to generate the matrix $\mathbf J^{-1/2}$, just take the eigenvectors from $\mathbf J$ (denoted $\mathbf Q$) and the inverse square root of the eigenvalues (denoted $\mathbf D^{-1/2}$) .. this is a "scalar"/"entrywise" operation because $\mathbf ...


6

Yes, for an SPD matrix $\mathbf A$ there are a variety of Cholesky-like decompositions, you can derive the $\mathbf A = \mathbf L^T \mathbf L$ variant by first writing down an educated/structured guess.. $\begin{bmatrix} \mathbf A_{11} & \mathbf a_{21}^T \\ \mathbf a_{21} & \alpha_{22} \end{bmatrix} = \begin{bmatrix} \mathbf L_{11}^T & \...


6

I don't think there is a way to display the used LAPACK function names natively during runtime using the interfaces provided by scipy.linalg. Depending on your goals you can: read the source code and deduce the logic from there. Unfortunately, this is not a runtime-use scenario, but a human analysis. fork your own version of scipy, add custom outputs (to ...


5

A related task to this is to find a subset of column vectors that are maximally linearly independent. Linear independence isn't exactly the same thing as asking for a large determinant, but if we can use one as a proxy for the other then this heuristic may help you. Rank revealing QR factorization (RRQR) is a good way to achieve this related task. One of ...


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 ...


4

It's NP-hard, but there are some heuristic algorithms. See https://arxiv.org/abs/1502.07838 which treats this problem.


4

I eventually found an answer. The image of a hypercube in $\Bbb R^3$ under a linear map is called a zonohedron. They can be calculated efficiently, for example using the algorithm in An Efficient Algorithm for Generating Zonohedra (PostScript file) by Paul Heckbert.


4

Since a SPD matrix is invertible, we can make the Cholesky decomposition $A^{-1} = PP^T$. Since $A$ is non-singular, so is $P$, and the inverse of a triangular matrix is triangular, so writing $L = P^{-1}$ we have $A^{-1}$ = $L^{-1}L^{-T}$. Inverting both sides gives $A = L^TL$.


4

You might want to use the fact that: $$ ||A||_2=\sigma_\max(A) $$ where $\sigma_\max$ is the largest singular value. If you are interested in details, this Math SO question should be interesting. Thus, $$ ||A^{-1}||_2=\frac{1}{\sigma_\min(A)} $$ where $\sigma_\min$ is the smallest singular value. You certainly want to avoid the actual calculation of the ...


4

I think you already know the answer- find the inverse of each of the small diagonal blocks. The details of how you might break the matrix up into the block in Python are something out of scope for this stackexchange group. As others have suggested you most likely don't really need the inverse of this matrix and could probably do everything you need with a ...


3

You don't need the inverse even with the goal of finding $K^{-1} h h^{T} K^{-1} - K^{-1}$. If you are interested to have this expression, I would explain how you can convert it to a matrix equation and then solve it more efficiently: Let's define the $X$ as: $$X = K^{-1} h h^{T} K^{-1} - K^{-1}$$ Your objective is to calculate $X$ in this equation by ...


3

I believe that $\mathbf W\cdot \mathbf 1$ denotes a matrix-vector product of a matrix $\mathbf W$ with a vector of ones $\mathbf 1$ of the corresponding size (number of columns). The result of this operation is a vector (of the size of number of rows of $\mathbf W$), so, it makes total sense to apply a diagonal function on it to arrive to a matrix form. ...


3

So, a cofactor matrix is a transpose of an adjugate matrix. I know of the following paper: G. W. Stewart, "On the adjugate matrix," Lin. Alg. Appl., vol. 283, no. 1–3, pp. 151–164, Nov. 1998. (available with full text) There, the author works on an algorithm of computing a adjugate matrix $\text{adj}(A)$ when $A$ is nearly singular or singular. For such ...


3

If it is known all the eigen values are positive, than a numerical brute force iterative method of the form $x_{t+1} = (1-c\mathbf{A})x_t$ where $c$ equals $1/\lambda_0$ where $\lambda_0$ is the dominant eigenvalue of $\mathbf{A}$ can be used. Expressing this more generally $x_{t+n} = (1-c\mathbf{A})^nx_t$ where $n$ represents the $n^{th}$ iteration. ...


3

Assuming that multiplication by a diagonal matrix is computationally easy, we can write it as: $$ ADA^T = (A\sqrt{D})(\sqrt{D}^TA^T) = BB^T $$ where $B=A\sqrt{D}$. Motivated by the fact that the left-singular vectors of $B$ are a set of orthonormal eigenvectors of $BB^T$, we can further proceed: $$ B = U\Sigma V^T $$ The first columns $m$ columns of $U$ (...


3

One can start by writing the eigenvalue problem $$ \left(A - \lambda_i\,I\right) w_i = 0, \tag{1} $$ with $\lambda_i$ one of the two eigenvalues and $w_i$ its eigenvector. By using the definition of $A$ $(1)$ can also be written as $$ \langle u, w_i\rangle u + \langle v, w_i\rangle v = \lambda_i\,w_i. \tag{2} $$ This implies that $w_i$ has to be a linear ...


3

The paper you linked to describes an algorithm similar to the Kabsch algorithm from what I see. It's used to find the least squares rotation between two sets of points. For your case you need something else entirely. Suppose sensor 1 has matrix $M_1$ at time t and sensor 2 has matrix $M_2$ at time t. That means that a vector, whose coordinates in the local ...


2

I believe that there is a good alternative that benefits from the geometry of the parameter space and completely eliminates the need for constrained optimization. If you explicitly wanted to make use of Lagrangians, I will definitely not be answering the question, but I thought it might be worthwhile to consider the perspective I will describe. In particular,...


2

Some thoughts on a particular case of low-rank matrices ($k=\text{rank}(M)\ll n)$. Here, I can suggest some economical version of SVDs: Rand-SVD Here is the R package documentation that also contains references for canonical papers of V. Rokhlin and P. Martinsson on this topic. That should reduce the complexity to $\mathcal O(nk^2+k^3)$ as opposed to $\...


2

You don't have to store $\tau$, you can re-compute it from the rest of the vector. (You can recompute $v_1$ from the other entries also in the normalized version, but it's clearly an unstable computation because of those subtractions.) Actually, you can reuse the lower-triangular part of $R$ to store $v_2,...v_n$, so that the factorization is computed fully ...


2

I don't think Octave's syms package supports this kind of matrix algebra operations, but you can do that in Mathematica. EDIT: and also in Sympy.


2

I was able to do it this way though it may had been easier to plot it the way Federico has mentioned before: imagesc(A). I declared an A matrix of zeros of mxn dimension and then loaded the ratings. I then applied HeatMap(A) which yields the image below. A = importdata('u.data'); user_id = A(:, 1); movie_id = A(:, 2); rating = A(:, 3); % Build matrix R ...


2

For Hadamard product, you have (source): $$ \text{rank}(A\odot X) \leq \text{rank}(A) \text{rank}(X) $$ This is as tight as you can go even with knowing $\text{rank}(X)$. Without knowing anything about $X$, there is certainly nothing you can say. (See the examples that nicoguaro and Wolfgang Bangerth brought in the comments section). Now, imagine you know ...


Only top voted, non community-wiki answers of a minimum length are eligible