Hot answers tagged

28

Stability under perturbations Let $E$ be a perturbation such that $\|E\| \leq \varepsilon$. If $A$ is symmetric, then the eigenvalues of $A+E$ are at a distance $\varepsilon$ from those of $A$. (Bauer-Fike Theorem.) If $A$ is non-symmetric, then the eigenvalues of $A+E$ could be just anywhere in the complex plane. No bound can be formulated a priori. (...


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


8

Amazing question which has a long answer, but I will try to be concise. In the context of Krylov subspace methods for general matrices, the eigenvalues of a non-symmetric matrix mean very little. In “Any nonincreasing convergence curve is possible for GMRES”, Greenbaum et al. show that any nonincreasing convergence curve is possible for GMRES independent of ...


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


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


6

Let's express the matrix $A \in \mathbb{R}^{n \times n}$ with which we want to solve linear systems as $$ A = S + U V $$ where $S$ is a symmetric matrix, $U \in \mathbb{R}^{n \times r}$, and $V \in \mathbb{R}^{r \times n}$. That is, $U V$ is a low rank update to account for the lack of symmetry. From your question, it appears $r$ is just 1 or 2. The ...


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

Usually, not always but usually, even if both $A$ and $x$ are sparse, $Ax$ is not. Even when it is, it is denser than $x$. If you consider something like $e^Ax$, which can be rewritten as $(I+A+A^2/2+\dots)x$, there is no guarantee that it will have the same sparsity pattern or the same number of nonzeros as $x$. Which means that either scipy would have to ...


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

Instead of using 4 levels of nested loops, you can take advantage of Kronecker products to simply your commutator_matrix function to def commutator_matrix(X): id = np.identity(np.shape(X)[0]) return np.kron(np.transpose(X), id) - np.kron(id, X) This is still the transpose of what you want. The function you are looking for is def ...


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

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


4

Another approach, which might be of interest to you is randomized sampling. This is of particular interest if you can quickly compute matrix-vector products $x\rightarrow Ax$ and $x\rightarrow A^* x$. The core idea is to form a small sampling matrix $S = A\Omega$, where $\Omega$ is a Gaussian random matrix. If the sampling matrix is large enough, $S$ will ...


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


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

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

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

Comment on: "But why are we so confident that we usually don't need to find the eigenvalues of non-symmetric matrix?" Eigenvalues and eigenspaces are relevant when you have an endomorphism, the same vector space as domain and range, like for instance in the linearization of an ODE system around some (equilibrium or periodic) solution. Then you can ...


3

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


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

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

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

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.


2

output[i*n + j][k*n + l] = com[k][l] That's your mistake I think -- reversed indices. To compute the matrix $M$ associated to a linear operator $f$ (the way it's usually taught in a linear algebra course), you need to take a basis of the input space $e_1, \dots, e_n$, compute $f(e_J)$ for each $J$, and write its coordinates (wrt a basis of the output space) ...


1

Did you try with Dask ? https://examples.dask.org/machine-learning/svd.html You can manage very large matrices. There is also a nice blog post about it https://blog.dask.org/2020/05/13/large-svds


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