66

The singular value decomposition, when applied to a real symmetric matrix $A = \sum_i \lambda_i(A) u_i(A) u_i(A)^T$, computes a stable mathematical object (spectral measure $\mu_A = \sum_i \delta_{\lambda_i(A)} u_i(A) u_i(A)^T$, which is a projection-valued measure) using a partially unstable coordinate system (the eigenvalues $\lambda_i(A)$ and eigenvectors ...

22

Recently there is a PhD thesis about the practical fast matrix multiplication algorithms like Strassen:
Matrix multiplication is a core building block for numerous scientific computing and, more recently, machine learning applications. Strassen's algorithm, the original Fast Matrix Multiplication (FMM) algorithm, has long fascinated computer scientists due ...

22

The SVD decomposition falls under the family of phenomena where discontinuity implies non-computability. (Intuitively, this is because at the point of discontinuity infinite precissions is required.)
In this particular case we speak of the (dis)continuity of a multivalued function which takes a matrix to any of its decompositions, or better, the non-...

18

This is primarily an issue of backwards vs. forwards stability. Good SVD algorithms are backwards stable in the sense that the computed singular values and singular vectors are the true singular values and singular vectors of a slightly perturbed problem. You may see this by noting that while $P$ may change drastically as you change $\epsilon$, the product $...

13

Here are a few relevant references (a precursor to "1." is also cited by Piyush Grover in the comments to the question):
Fast low-rank modifications of the thin singular value decomposition by M. Brand, Linear Algebra and its Applications, 415 (2006)
Restricted rank modification of the symmetric eigenvalue problem: Theoretical considerations, P. Arbenz, W. ...

13

Implementing the QR factorization with Householder rotations is cheaper ($2n^2m$ vs $3n^2m$ for a $m\times n$ matrix), and equally accurate in practice. See Section 19.6 of Higham's Accuracy and Stability of Numerical Algorithms, or Golub-Van Loan for more explicit algorithms.
Moreover, in a Householder-based implementation there is a higher fraction of ...

13

Disclaimer 1: Treating these topics properly would require a quick course in numerical analysis.
Disclaimer 2: If you are using any sane computer system, it's already going to have a library function to solve linear systems implemented, which is going to be better of what you can code yourself if you don't have a solid grasp of numerical linear algebra and ...

answered May 11 '19 at 15:14

Federico Poloni

15.6k22 gold badges6161 silver badges9898 bronze badges

12

@ user40484 , fortunately your estimate for the complexity of SVD is not optimal. Otherwise, you put unemployed specialists in image compression. The complexity is in $O(\min(mn^2,m^2n))$.
Assume the data points are in the columns of $A\in M_{m,n}(\mathbb{R})$ where $m\leq n$. Note that $AA^T$ is the dataset covariance matrix. Then a simple method is to ...

12

The exponential of a matrix is not, in general, computable explicitly. Strictly speaking, the Complexity of Matrix Exponential does not exist. Instead, you have an approximate calculation and you are asking about its complexity. Now, a few comments :
there is not just one approximate exponential, but actually a sequence of approximations. Therefore, the ...

answered May 17 '16 at 14:27

Denis Serre

40.3k88 gold badges109109 silver badges225225 bronze badges

11

Not in general. An explicit and elementary counterexample
is the sparse triangular matrix with $1$'s on the diagonal
and $-1$'s just above it: the inverse is the triangular matrix with
every entry on or above the diagonal equal $1$.

matrices na.numerical-analysis approximation-algorithms nonlinear-optimization numerical-linear-algebra

answered Aug 19 '13 at 17:37

Noam D. Elkies

63k1212 gold badges220220 silver badges299299 bronze badges

10

Moler's paper "Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later" contains the following extracts:
In estimating the time required by matrix computations it is traditional to estimate the number of multiplications and then employ some factor to account for the other operations. We suggest making this slightly more ...

10

There are polynomial-time reductions from your problem to Graph Isomorphism and vice-versa.
As a quick definition, when I speak of 'subdividing' an edge, I mean to replace each edge $u, v$ with a path $u, w, v$ where $w$ is a new 'midpoint' vertex.
Transposable $\rightarrow$ GI:
Convert your $n \times n$ matrix into a bipartite graph $H$ with coloured ...

answered Nov 14 '18 at 18:27

Adam P. Goucher

10.3k22 gold badges4242 silver badges8080 bronze badges

9

The problem is NP-hard over any field not of characteristic 2. We show this by reduction from an NP-complete NAE-3-SAT problem of checking satisfiability of a boolean formula $\land_{i} \mathrm{NAE}(z_{i, 1}, z_{i, 2}, z_{i, 3})$, where $z_{i, j}$ are variables $x_{\cdot}$ or their negations, and $\mathrm{NAE}(\ldots)$ is true iff not all of its arguments ...

9

For an $n\times n$ matrix, the answer is
$$|\det A|^{\frac1n}.$$
Explanation: on the one hand, $\rho(UA)\ge|\det (UA)|^{\frac1n}=|\det A|^{\frac1n}$. On the other hand, singular value decomposition gives $A=PDQ$ where $P,Q$ are orthogonal and $D>0$ is diagonal. Then
$\rho(UA)=\rho(QUPD)$ and
$$\min_U\rho(UA)=\min_V\rho(VD)$$
where $V$ runs over the ...

8

I think Pushpendre's answer isn't quite right, but it gets you most of the way there. Getting rid of that pesky constant term is a bit tricky relative to the homogeneous case.
Let's take his suggested substitutions:
$$
\begin{array}{rl}
M&:=N^\top N\\
y&:=Nx\\
D&:=CN^{-1}\\
B&:=N^{-\top}AN^{-1}
\end{array}
$$
I don't think these are ...

8

The linear-in-$g$ system is uniquely solvable whenever $f_0\ne0$. therefore the solution set is parametrized by $f$ in the form
$$g=f_0^{-10}P(f;a),$$
where $P$ is linear in $a$, polynomial in $f$ and homogenenous of degree $9$. Not only it is a continuum, but it has dimension $10$.
When $f_0=0$, there is no solution in general, say for $a_1\ne0$.

7

Just wanted to note that if one knows that the rank two update is of the form $uv^t + vu^t$ and one knows $v$ and $u$ then it is easy to find $x,y$ such that $uv^t + vu^t = xx^t - yy^t$.
$x,y = \sqrt{\frac{|v|}{2|u|}} (u \pm \frac{|u|}{|v|}v)$. Once we know $x,y$ then the problem reduces to two symmetric rank one updates and we can use the symmetric rank ...

7

That matrix is the product of two positive definite matrices, $M=UW$, with $W=A^TVA$. Hence it is similar to the symmetric matrix $U^{1/2}WU^{1/2}$, which is congruent to $W$ and thus positive definite. For a slightly more general result, see Horn, Johnson, Matrix analysis, 1st ed, Theorem 7.6.3.

7

For a short fat matrix $G$ (more columns than rows), $\|AG\|_F \geq \sigma_{\min}(G)\|A\|_F \geq n \sigma_{\min}(G) \|A\|$, where $\sigma_{\min}(G)$ is the least singular value of $G$. This follows from the fact that $\|G^*v\|\geq \sigma_{\min}(G^*)\|v\| = \sigma_{\min}(G)\|v\|$ for every vector $v$, applied to the rows of $A$.
However, if $G$ is tall thin (...

answered Nov 17 '16 at 13:06

Federico Poloni

15.6k22 gold badges6161 silver badges9898 bronze badges

7

To make sure I understand: you have an $m \times n$ matrix and vector in $\mathbb{Z}^n.$ Is the system a priori over- or under-determined? The former case is a little easier than the latter, but in any case, your problem is the so-called intermediate expression swell, and the easiest way to get around it is to solve it mod many primes and then chinese-...

7

With the traditional algorithms and complexity measures used in numerical linear algebra (dense real matrices, floating point computations, flop count as a complexity measure), they are both more or less equally fast, with the characteristic polynomial probably being slightly faster.
Eigenvalue computation requires an iteration so its time is not constant, ...

answered Jul 19 '19 at 21:47

Federico Poloni

15.6k22 gold badges6161 silver badges9898 bronze badges

6

I came back to this a few months ago and I can now answer my own question. I hope it is appropriate to answer my own question given the length of time.
Bernhard Beckermann and I just submitted a paper [1] that shows that if $AX-XB = F$ with $A$ and $B$ normal matrices, then the singular values of $X$ can be bounded as
$$
\sigma_{1+\nu k}(X) \leq Z_k(\sigma(...

6

While not a full characterization, the following result on $\sigma_n$ shows that $\sigma_k > \epsilon_n$, since $\sigma_n \ge \epsilon_n$.
Following my answer on this MO question here, we see that
\begin{equation*}
\sigma_n \ge \frac{\det(H)}{\left[m + s/\sqrt{n-1}\right]^{n-1}} =: \epsilon_n,
\end{equation*}
where
\begin{eqnarray*}
m &=& \...

6

To answer the question: "...is there a name for this matrix..."
Beyond the "Gaussian-Toeplitz matrix" mentioned in the comments, the said matrix is a special case of the Gaussian Kernel (which is also the alluded to "heat kernel"):
\begin{equation*}
k(x_i,x_j) = e^{-c(x_i-x_j)^2}\qquad x_1,\ldots,x_n \in \mathbb{R}.
\end{equation*}
Clearly, the matrix $A =...

6

Your problem has been answered at
https://scicomp.stackexchange.com/questions/14096/sparse-smallest-eigenvalue-problem-on-a-linear-subspace :) Or you can read Golub's original paper Some modified matrix eigenvalue problems
The basic intuition is that basically you want to find the eigenvalues over a subspace that is defined $Cx = b$ so just find its basis ...

6

The set of $n \times n$ symmetric Toeplitz matrices is
$$\left\{ x_1 \mathrm M_1 + x_2 \mathrm M_2 + \cdots + x_n \mathrm M_n \mid x_1, x_2, \dots, x_n \in \mathbb R \right\}$$
where $\mathrm M_1, \mathrm M_2, \dots, \mathrm M_n$ are $n \times n$ symmetric Toeplitz basis matrices. Let $\mathrm M_1 = \mathrm I_n$ correspond to the main diagonal, whereas the ...

answered Mar 14 '18 at 21:39

Rodrigo de Azevedo

2,02233 gold badges1111 silver badges2525 bronze badges

6

Removing all unnecessary parameters, we come to the equation $\Omega^{-1}=2 W\odot \Omega + B$ where $B$ is positive definite. We need to find a solution in the cone $M_+$ of positive definite matrices. The solutions are exactly the stationary points of
$F(\Omega)=\log\det\Omega- \operatorname{Tr} [(W\odot\Omega)\Omega+B\Omega]$ (I hope that is not where ...

6

Clearly, a necessary condition is that for every word $w$ in two letters, one has
$${\rm Tr}\,w(A^t,B^t)={\rm Tr}\,w(A,B).$$
Equivalently,
$${\rm Tr}\,\hat w(A,B)={\rm Tr}\,w(A,B),$$
where $\hat w$ is the reverse word. Namely, if $w=x^\ell y^mx^n\cdots$, then $\hat w=\cdots x^ny^mx^\ell$.
Unless the word $\hat w$ be conjugated (in the free group ${\mathbb F}...

answered Nov 15 '16 at 10:48

Denis Serre

40.3k88 gold badges109109 silver badges225225 bronze badges

5

You can use the Chebyshev Polynomial expansion
to calculate the effect of the matrix exponential on a vector.
Which is a standard technique in quantum chemistry community and the method is extremely stable and fast. This method was developed by Tal-Ezer and Kossloff in an article named An accurate and efficient scheme for propagating the time dependent ...

5

Just a followup to @Lucia's comment, here is the paper:

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

#### Related Tags

numerical-linear-algebra × 204linear-algebra × 128

matrices × 68

na.numerical-analysis × 65

eigenvalues × 19

matrix-analysis × 16

matrix-theory × 10

sparse-matrices × 10

algorithms × 9

nonlinear-optimization × 9

matrix-equations × 8

fa.functional-analysis × 7

oc.optimization-and-control × 7

convex-optimization × 7

linear-programming × 7

eigenvector × 7

tensor × 6

reference-request × 5

multilinear-algebra × 5

matrix-inverse × 5

computational-complexity × 4

random-matrices × 4

perturbation-theory × 4

ra.rings-and-algebras × 3

tensor-products × 3