Hot answers tagged

12

If you can compute products with $A$ and $A^T$, as you specify in a comment, you can run the classical sparse SVD algorithms such as scipy.sparse.linalg.svds, Matlab's svds, or Julia's Arpack.svds, which are based on Lanczos bidiagonalization. They are designed to compute singular values, and are likely to be more robust than a minimization routine coded by ...


9

An anti-Hermitian matrix is diagonalizable, with orthogonal eigenvectors (ref). Hence you can write $X = PDP^{-1}$, where $D$ is a diagonal matrix. Therefore the exponential can be calculated as $e^X=Pe^DP^{-1}$, and $e^{tX} = Pe^{tD}P^{-1}$. If $d_1$, $d_2$, etc.... are the diagonal elements of $D$, then for each value of $t_i \in t$, $e^{t_iD}$ is just a ...


7

This is a slightly modified version of my response on math.stackexchange. One standard approach to computing matrix functions times a vector $f(M)x$ or quadratic forms $x^Tf(M)x$ when $M$ is symmetric is via the Lanczos algorithm. Lanczos computes an orthonormal basis $Q_k = [q_1, \ldots, q_k]$ for Krylov subspace $\operatorname{span}(x,Ax,\ldots, A^{k-1}x)$ ...


7

This is my comment expanded into an answer. MATLAB is a column-major programming language; it is not very hard to find it in the documentation once you know what to look for: https://www.mathworks.com/help/matlab/matlab_external/matlab-data.html#f22019 The layout of data structures (most commonly row-major or column-major) used in a programming language is ...


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


4

Before even discussing any vectorization, there are gross inefficiencies in the code you provided. Just by employing basic techniques, I got 21-fold speedup: n = 1000000; y = zeros(1,n); tProbs = [0.1 0.5; 0.9 0.5]; iProbs = sum(tProbs^100000,2); iProbsSqr = iProbs.^2; V0 = [1;1]; for i=1:n V0 = tProbs*V0; y(i) = sqrt(dot(V0.^2,iProbsSqr)); end If the ...


4

tch already gave you a good answer, but I have a suggestion for a simpler problem that you could address pretty easily and which should help you sanity check your results. Let $$M = \sum_i\lambda_iv_iv_i^\top$$ be the eigenvalue decomposition of $M$ and furthermore take the eigenvalues to be sorted in decreasing order: $\lambda_1 \ge \ldots \ge \lambda_n$. ...


4

LAPACK has an implementation of the svd of a 2x2 triangular matrix. It appears to be very robust. The routine is XLASV2. To apply to a regular 2x2 matrix, you can simply apply a single givens rotation from the left/right.


4

For symmetric, more generally, normal, matrices there is no danger in diagonalizing a matrix to evaluate a matrix function, because all eigenvectors are perfectly conditioned. You should get answers with errors close to the machine epsilon times the condition number of the function evaluation problem. This is as accurate as you can get in floating point. ...


4

I have not quite understood what the dimensions of your matrices are (see comments). but the following code may give you a start. It is not efficient or particularly elegant, but you can improve from there: import numpy as np A = np.matrix('1 2; 3 4') B = np.matrix('1 2; 3 4') C = np.matrix('1 2; 3 4') result = (np.matmul(A.transpose(),np....


3

Yes, the eigenvectors found with this method may depend on $x$ and $y$, but no, it doesn't matter in practice. If $A$ and $B$ share a basis of common eigenvectors, then $$ A = V\operatorname{diag}(\alpha_1,\dots,\alpha_n)V^{-1}, \quad B = V\operatorname{diag}(\beta_1,\dots,\beta_n)V^{-1}. $$ If $A$ has all distinct eigenvalues, then $V$ is unique, up to ...


3

This problem is known as joint diagonalization, and it has two variants: orthogonal, in which the basis vectors are orthonormal, and non-orthogonal, which is harder to solve, but which may be more appropriate to your application. The simplest method I know of seeks a unitary matrix $U$ that minimizes the sum of squares of the off-diagonal elements of $U^HAU$...


3

If your advection problem had Dirichlet of Neumann boundary conditions, the linear system would be tridiagonal and you could apply the Thomas algorithm. With periodic boundary conditions, however, we lose this. If c(x) is a constant independent of x, the matrix would be circulant and linear systems could be solved efficiently using FFTs. An even better ...


3

The Frobenius norm of a binary matrix is the square root of the number of non-zero elements. Let the point (0.,...0) be the origin, and let's say the vec'd binary matrix elements are the coordinates of a point, with coordinates being zero or one. Then the Frobenius norm is the Euclidean distance from the origin to that point.


3

Another nice way to do this sort of thing is numpy.einsum, which allows you to express the multiplication in index notation: >>> A=np.array([[1,2,3],[4,5,6]]) >>> B=np.array([[1,2],[3,4]]) >>> C=np.array([[1,2,3],[4,5,6]]) >>> A.T@B@C array([[ 85, 116, 147], [113, 154, 195], [141, 192, 243]]) >>> np....


3

mats_row=np.array([[A1,A2,A3,...,A_n]]) #Create array of matrices with shape (1,M,N,N) mats_column=np.transpose(mats_row,(1,0,2,3)) #Make a copy with shape (M,1,N,N) block=.5*(mats_row+mats_column) #Add, broadcasting to a (M,M,N,N) array This works by utilizing Numpy's broadcasting capabilities. The basic idea is similar to if you added a matrix where each ...


3

A dyadic product takes as input two vectors and outputs a second order tensor. This is what I know as a dyadic product, and a dyad is the term $\mathbf{a}\mathbf{b}$. A general second order tensor can be written as a linear combination of dyads. Commonly the symbol $\otimes$ is referred as the tensor product and it outputs higher-order tensors. I think that ...


2

I am not sure about your application -- and we say the $L^2$ norm of a function and not a system. But for simplicity I will explain the concepts for real valued functions. Consider an open domain $\Omega$ and a function $f:\Omega \to \mathbb{R}$. We say that $f \in L^2(\Omega)$ if $||f||_{L^2(\Omega)} < \infty$ where \begin{equation} ||f||^2_{L^2(\Omega)} ...


2

The first equality you wrote is not correct, as noted by other users. However, what I think you want to know is why $(I:A)I = (I \otimes I)A$. You can show this just by using dyads properties. $$(B \otimes B) : A =(B_{ij}B_{kl} e_i \otimes e_j \otimes e_k \otimes e_l):(A_{mn} e_m \otimes e_n) = B_{ij}B_{kl} A_{mn} \delta_{km} \delta_{ln} e_i \otimes e_j = B_{...


2

This landmark paper by George proves that a nested dissection ordering of a regular, four-node element, finite element mesh produces minimum fill-in. Although it is straightforward to produce such an ordering by inspection, the graph algorithms that attempt to do this for a general sparse structure only approximate this ordering. Assuming you are obtaining ...


2

The inverse of the (1,1) block of $$ \begin{bmatrix} A & B\\ C & D \end{bmatrix}^{-1} $$ is $A-BD^{-1}C$ (Schur complement). This is what you are trying to compute, if I understand correctly from your explanation ("marginalize" may be standard in your domain, but it is not standard linear algebra language). So at least you can reduce to ...


2

@Federico Poloni 's fine answer states the impossibility of getting an exact yes/no answer using IEEE arithmetic. However, using interval arithmetic with outward rounding, it is possible to get a "not singular/don't know" answer. In particular, it may be possible to definitively conclude that the smallest singular value is strictly greater than ...


2

I also suggest looking into the condition number estimators, which will (with some degree of [un]reliability) predict how effectively numerically singular the matrix is. In particular, "Spectral Condition-Number Estimation of Large Sparse Matrices" based on LSQR seems like an interesting choice. attracted my attention one day. I would also ...


2

We can do some transformations to your problem to show that it's easily solvable via a linear program: Given a matrix $M$ with non-negative real entries and a vector $v$ you wish to solve the problem: $$ \begin{align} \min_v \quad & \lVert Mv \rVert_\infty \\ s.t. \quad & v_i\ge0 \\ & \sum_i v_i = 1 \end{align} $$ Now, note that $\lVert Mv \...


2

$\def\p{\partial_{p}}$The linked paper has not been peer-reviewed and the result is clearly wrong. Expanding the expression in index notation yields $$\eqalign{ {\cal H} &= \nabla(A\cdot B) \\ {\cal H}_{pik} &= \p(A_{ij}B_{jk}) \\ &= (\p A_{ij})B_{jk} + A_{ij}(\p B_{jk}) \\ &= (\p A_{ij})B_{jk} + (\p B_{kj}^T)A_{ji}^T \\ {\cal H} &= (\...


2

As its documentation suggests in multiple places, rref is "mainly of academic interest" (read: "used only to explain Gaussian elimination to undergrads"), and is not a serious competitor of SVD-based algorithms in terms of stability. I recommend against it.


1

Probably you accumulated rotations wrong. At each step, you have a working matrix $D$ that starts from $D=A$ and should converge to a diagonal matrix, and a matrix $Q$ that accumulates all rotations performed on $D$. Print after each step $\|QDQ^T - A\|$, and see where it stops being small.


1

Here's my 2 cents. I would set up the following minimization problem $$ \pi(x) = \frac{1}{2} (Ax)^T(Ax) $$ If $A$ has eigenvalues which are zero, there will exist a nonzero $x$ such that $\pi(x)=0$. So, I would try computing the gradient of $\pi$ wrt $x$ and use a gradient descent algorithm to drive $\pi$ towards zero. If you get reasonably close to zero ($\...


1

If you're ok with positive semi-definite matrices, then you can refer to this paper, which I think is what Frederico's answer also explains: Higham NJ. Computing a nearest symmetric positive semidefinite matrix. Linear Algebra and its Applications. 1988 May;103(C):103-118. A python implementation can be found in this answer, which also provides a link to a ...


1

As mentioned, netlib BLAS is not at all optimized, but it is definetly the "refblas". Using IKML, ACML, OpenBLAS or "your vendor" BLAS, you are (somehow) assured, that the results of the operation of the optimized BLAS is equal to the "refblas" up to a known error. Take into care that: vendors (intel, amd, nvidia, ...) try hard ...


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