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


6

Let me see if I can restate the problem to answer your question: You have 5 sets of data: $y_{i}$, $x_{i}$ such that for each $i = 1, \ldots, 5$, $x_{i}$ and $y_{i}$ are both 3 by 1 vectors. You also have a 3 by 3 matrix $B$ that is supposed to relate $x_{i}$ to $y_{i}$ via: $$y_{i} = Bx_{i} $$ Is this formulation correct? If so, the proper methods for ...


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

SLICOT is the tool to use for dense problems. For large but sparse system, there is the lyapack toolbox for MATLAB. The algorithms in lyapack base on computing iteratively low-rank factors $Z_n$, so that $Z_n^HZ_n$ approaches $\Sigma$, where $\Sigma$ is the symmetric (positive or negative) definite solution of the Lyapunov equation. Computing only the ...


5

One of the best-performing methods for solving large-scale Lyapunov equations is ADI. It is an iterative algorithm that returns an approximate low-rank decomposition $X \approx VV^T$ of the solution $X$. In this case, you can work with this decomposition of both Gramians to reduce the eigenproblem to a smaller one. I suggest you to start approaching this ...


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

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

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

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


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

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

Generalizing Arnold's example, suppose $D = I$. Then the problem is to find an orthogonal $X$ satisfying $X(a-b) = 0$, which is possible only when $a = b$. Interestingly, in 2D with $D = \text{diag}(1, -1)$, one can always find a plane rotation $X$ that does the job. Evidently the solvability for arbitrary unit-norm $a$, $b$ depends only on $D$. What is ...


3

You could connect to MATLAB using this. Your matrices aren't too large : hand coding the algorithms shouldn't result in too much time loss, maybe it'll run for 1 hour. It may or may not be too long depending on various factors. Though, coding it yourself may not be easy at all. I don't think I can, and I've been dealing with this for the past few months. ...


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

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

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


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

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


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

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


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


1

Your belief that there should always be a solution is wrong for $n=1$ if $D=a=1$ and $b=-1$, which satisfies all your requirements. Thus it seems unlikely that one can say anything in general. Posing the linear and the orthogonality constraint as a least squares problems in the Frobenius norm, and submitting it to an unconstrained optimization routine is ...


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