8

It is called a T-Sylvester equation, or *-Sylvester equation in the complex case. Solvability conditions and a pseudocode algorithm based on the Schur form are in https://doi.org/10.13001/1081-3810.1479 . Analogous considerations for a more general class of equations and a Fortran-90 implementation of the last step of the resulting solution algorithm (the ...


7

As was mentioned in the comment, calculating $x=M^{-1}y$ is equivalent to solving $Mx=y$. Here is the full solution: First, you can reformulate the equation to: $Bx=(2A+I)(C^{-1}+A)b$, and by defining $\tilde{b}=C^{-1}b$, the equation can be rewritten as: $Bx=(2A+I)(I+AC)\tilde{b}$. First, compute $\tilde{b}$ by solving (using, for example, LU ...


7

Took me quite a while to figure this out and as usual it becomes obvious after you find the culprit. After checking the problematic cases reported in David S. Watkins. A case where balancing is harmful. Electron. Trans. Numer. Anal, 23:1–4, 2006 and also the discussion here (both being cited in arXiv:1401.5766v1), it turns out that matlab uses the ...


5

$$\rm x = B^{-1} (2A + I) (C^{-1} + A) b$$ Left-multiplying both sides by $\rm B$, $$\rm B x = (2A + I) (C^{-1} + A) b$$ Let $\rm y:= C^{-1} b$. Hence, $$\rm B x = (2A + I) (C^{-1} b + A b) = (2A + I) (y + A b) = (2A + I) y + (2A + I) A b$$ and, thus, we obtain a linear system of $2n$ equations in $2n$ unknowns $$\begin{bmatrix} \mathrm B & -(2 \...


5

I will try to give my thought on the first question regarding fast $3\times 3$ inverse. Consider $$ A=\left[ \begin{array}{ccc} a & d & g\\ b & e & h\\ c & f & i \end{array}\right] $$ Since the matrices are small and very general (do not feature any known structure, zeroes, relative scales of the elements), I think it would be ...


5

I recommend ARPACK which provides matrix-free routines for generalized eigenvalue problems with examples in their documentation for this purpose. This method is one of the most widely used, especially for relatively large problems where you are only searching for a few eigenvalues/vectors (such as the smallest, in your case). Further, eigenvalue problems ...


5

There are specialized methods for the minimization of a differentiable function $f(X)$ subject to the orthogonality constraint $X^{T}X=I$. See for example: Lai, Rongjie, and Stanley Osher. “A Splitting Method for Orthogonality Constrained Problems.” Journal of Scientific Computing 58, no. 2 (2014): 431–449. Wen, Zaiwen, and Wotao Yin. “A Feasible Method ...


4

The beauty of iterative methods is that all they require you to do is matrix-vector multiplications. In your case, the product of your matrix $A$ with a vector $y$ can be written as $z=Ay = (B^TCB)y= B^T(C(By))$ which shows that all you need is three matrix-vector products but no matrix-matrix products. If $C$ happens to be symmetric and positive semi-...


4

$$\big\| \mathrm A \begin{bmatrix} \mathrm X & \mathrm X^2\end{bmatrix} - \begin{bmatrix} \mathrm B_1 & \mathrm B_2\end{bmatrix} \big\|_{\text{F}}^2 = \underbrace{\| \mathrm A \mathrm X - \mathrm B_1 \|_{\text{F}}^2}_{=: f_1 (\mathrm X)} + \underbrace{\| \mathrm A \mathrm X^2 - \mathrm B_2 \|_{\text{F}}^2}_{=: f_2 (\mathrm X)}$$ Everybody knows that ...


3

A simple way to introduce orthogonality constraints is to parametrize all orthogonal matrices using, either, the Cayley transform, $Q=(I-A)(I+A)^{-1}$, or the matrix exponential, $Q = \exp(A)$. In both cases, $Q$ will be orthogonal if $A$ is skew-symmetric. Searching over the space of skew-symmetric matrices is easy, as an element of the space can be ...


3

I encountered a similar problem in the past and I could find no simple solution either. One of the terms is a Kronecker product, another is a rank-1 modification, but the rest makes the problem more difficult. I don't think there is a closed-form solution; you could try using an iterative method, dropping some terms to get a preconditioner. But if someone ...


3

We have a system of two equations in $\mathrm X \in \mathbb R^{n \times n}$ $$\begin{aligned}\rm D X a &= \rm X b\\ \rm X^\top X &= \rm I_n\end{aligned}$$ The convex hull of the orthogonal group $\mathrm O (n)$ is defined by $\mathrm X^\top \mathrm X \preceq \mathrm I_n$, or, equivalently, by the inequality $\| \mathrm X \|_2 \leq 1$. Hence, we ...


3

SLICOT's algorithm is not that complicated, it's a reduction to Schur form + some back-substitution. You can check the Bartels-Stewart paper http://dl.acm.org/citation.cfm?id=361582 which is reasonably readable and explains how it works. The paper is about the nonsymmetric case, but it shouldn't be hard to adapt it to the symmetric one --- you just need one ...


3

The procedure you need to get the "reduced" equations is often referred to as "static condensation" in the FEM literature. You can partition your FE equations as follows: $$ \left[\begin{array}{cc} K_{cc} & K_{ci} \\ K_{ic} & K_{ii} \\ \end{array}\right] \left\{\begin{array}{c} V_c \\ V_i \end{array}\right\} = \left\{\begin{array}{c} Q_c \\ Q_i \...


2

Assume the equations are discretized on the $\tau$ grid, and introduce several column-vectors, using the notation where the superscript stands for the grid point index i $\in$ [1,...,n]. $\vec{\tau}=\left[ \tau^1,\tau^2,...,\tau^{n} \right] $ $\vec{\phi_1}=\left[ \phi_1^1,\phi_1^2,...,\phi_1^{n} \right] $ $\vec{\phi_2}=\left[ \phi_2^1,\phi_2^2,...,\phi_2^{...


2

Ignoring Newton's method here is the wrong approach! The fact that you're using Newton's method is what makes this cheap to add, and is what makes singular mass matrices possible. Essentially look at Implicit Euler. You'll notice that the only matrix you need to use in your quasi-Newton steps is $$ W = (I - \gamma J) $$ and if you have a mass matrix, this ...


2

Let us consider the simplest case, i.e., the case where $n = 1$. Rephrasing slightly: Given symmetric and positive definite matrix $\mathrm A \in \mathbb R^{m \times m}$ positive semidefinite diagonal matrix $\mathrm D \in \mathbb R^{p \times p}$ find a tall matrix $\mathrm X \in \mathbb R^{m \times p}$ with orthonormal columns such that $\rm A \approx X ...


2

If $N$ is on the order of 100,000 and $m$ is on the order of $100$, Then $J$ requires about 80 gigabytes to store in double precision and $V$ requires a trivial amount of storage. The product $M=JV$ is of size $N$ by $m$ and would be fully dense, requiring about 80 megabytes to store. You should have no trouble storing $V$ or $M$, but you may have to keep ...


2

You can prove convergence by satisfying the spectral radius relationship you note, choosing $S$ and $T$ such that $\rho(S^{-1}T) < 1$. This comes about by first writing two equations based on your operator splitting: $$ Sx = Tx + b$$ $$ Sx^{(k+1)} = Tx^{(k)} + b $$ where $x$ is exact solution and $x^{(k)}$ is the $k^{th}$ iteration's solution. Now ...


2

To expand on Bill Greene's answer, Fortran uses different syntax for single and double precision floating point numbers. 1.e0 is a single precision representation of the number '1.0', whereas 1.d0 is double precision representation of the number '1.0'. If you have a rather straight forward problem, and an associated low condition number, you will not notice ...


2

The matrix in your case 2 example has a much higher condition number than the one in the original 5-equation example. Accordingly, errors in the definition of the matrix terms result in larger errors in the solution. I believe that if you output the matrices from MATLAB with 16 digits of precision, the fortran codes will yield the same solution as MATLAB. ...


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


1

I wound up using a different approach to get the optimization working. In the paper, they use an approach based on Jacobi rotations where they rotate pairs of modes by a calculated optimal angle to increase the function value. This approach had previously been used to perform a similar orthogonal optimization to localize molecular orbitals (Rev. Mod. Phys. ...


1

You can always try a conjugate gradient (https://en.wikipedia.org/wiki/Conjugate_gradient_method). Once you have found the solution via LU, if the new changes affect only to a small number of equations (or even a large one), the convergence should be very fast. Here it is clearly explained how to do it: ftp://ftp.numerical.rl.ac.uk/pub/talks/isd_stanford50....


1

You are right that writing an efficient and robust eigenvalue solvers for large sparse systems is a difficult task that should only be done yourself as the very last resort. Spencer already gives the major players in his answer (ARPACK for single-node multithreaded computing, which is part of MATLAB and SciPy, and SLEPc for distributed computing, which can ...


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