Skip to main content
Share Your Experience: Take the 2024 Developer Survey
8 votes
Accepted

Solving $AX+X^TB=C$?

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....
Federico Poloni's user avatar
8 votes

Compute $x = B^{-1}(2A+I)(C^{-1}+A)b$ without calculating matrix inverses

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 ...
Gil's user avatar
  • 392
7 votes
Accepted

Matrix Balancing Algorithm

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 ...
percusse's user avatar
  • 393
5 votes
Accepted

Maximize a function of an orthogonal matrix

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 ...
Brian Borchers's user avatar
5 votes

Recurrence relation for matrices

This is an instance of the Riccati equation, which can be solved using the Matrix Sign function. Relevant section from Higham's "Functions of Matrices" book:
Yaroslav Bulatov's user avatar
5 votes

The fast, and The Backward-Stable (left) $3\times 3$ matrix inverse

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}\...
Anton Menshov's user avatar
  • 8,672
5 votes

Confusion about matrix differentiation in a nonlinear matrix equation

If a matrix is differentiated with respect to itself, the result should be a fourth order tensor. The easist way to see this is to work with components. $$ \frac{ \partial K_{ij}}{\partial {K_{kl}}} = ...
NNN's user avatar
  • 760
4 votes
Accepted

Mass Matrix and how to handle it (ODEs) - References

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 ...
Chris Rackauckas's user avatar
4 votes
Accepted

Solving linear system with matrix multiplication

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= ...
Wolfgang Bangerth's user avatar
4 votes
Accepted

Solving underdetermined Lyapunov equation?

Since $A$ is symmetric, it has an eigendecomposition $A = QDQ^*$ with $Q$ orthogonal. Then $$ M = A \otimes I + I\otimes A = (Q\otimes Q)(D\otimes I + I \otimes D)(Q\otimes Q)^* $$ is an ...
Federico Poloni's user avatar
3 votes
Accepted

Recurrence relation for matrices

Another different way to obtain an equivalent formula: let $Y = A^{1/2}XA^{1/2}$. Then, multiplying your equation from both sides by $A^{1/2}$, we have $Y^2+Y = A^{1/2}BA^{1/2} = C$. The matrix $C$ ...
Federico Poloni's user avatar
3 votes

Solution to Sylvester-like equation

EDIT: this answer is essentially useless after the question was updated. If the coefficients $A,B$ are symmetric, these approaches reduce to the simpler closed formula with the eigenvalue ...
Federico Poloni's user avatar
3 votes
Accepted

Solution to Sylvester-like equation

Working with a more standard form of a Stein equation, sometimes called (generalized?) Discrete Lyapunov Equation $$AXB - X = C$$ Rewrite it in standard linear equation form, use least squares solver,...
Yaroslav Bulatov's user avatar
3 votes
Accepted

Matrix regularisation for ill-conditioned problems

matrix regularisation can improve the stability of LU or Cholesky decomposition of ill conditioned problems. Not really, at least in the way the word "stability" is typically used in the ...
Federico Poloni's user avatar
3 votes

Solving MX=N where M is structured as a Gaussian 4th-moment tensor

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 ...
Federico Poloni's user avatar
3 votes

Maximize a function of an orthogonal matrix

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 ...
Amit Hochman's user avatar
  • 1,081
3 votes

Invert a matrix only on a subset of variables / Compute the "equivalent circuit"

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} ...
Bill Greene's user avatar
  • 6,094
3 votes
Accepted

Solution $X$ for $X(X^TX)^{-1}=X(Y^TY)^{-1}$

Notice that $X(X^T X)^{-1}$ is the Moore-Penrose inverse of $X^T$. Using properties of pseudo-inverse we get $(Y^TY)(X^T X)^\dagger = I$ For instance $X=Y$ or $X=\text{Cholesky}(Y^T Y)$
Yaroslav Bulatov's user avatar
3 votes

How can I express the solution to a discrete Lyapunov equation as an eigenvalue problem?

For completeness, I'm including a full derivation for the solution of the discrete Lyapunov equation that is adapted from @FedericoPoloni's paper (specifically section 4.1). For the derivation for the ...
mhdadk's user avatar
  • 165
3 votes

Matrix derivative

Observing that both the trace and the multiplication with a matrix $A$ are linear operations, it is easy to apply the chain rule. For this, first see that $$ \frac{\partial}{\partial X_{ij}} [\text{...
Wolfgang Bangerth's user avatar
3 votes
Accepted

Matrix derivative

Per the very helpful wikipedia page "derivative of the exponential map", based around the "Baker-Campbell-Hausdorff formula", the relevant formula for a directional derivative is $\...
Lutz Lehmann's user avatar
  • 6,109
3 votes

Solve linear system for only part of the solution vector

The only way that I'm aware of to take advantage of only needing a partial solution is in the triangular solves. But those are usually a small part of the total time in dense LU (unless you have many ...
Neil Lindquist's user avatar
2 votes
Accepted

Least Squares with Dense-Block Diagonal Structure

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$ ...
Brian Borchers's user avatar
2 votes
Accepted

Convergence conditions of a stationary iteration method for linear systems

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 ...
spektr's user avatar
  • 4,248
2 votes

Solve $AX=B$ where $A$ is a skyline matrix

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 ...
cbcoutinho's user avatar

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