Derivative of cholesky factorization

[blegat]

Given a Cholesky factorization Q = U'U and a symmetric matrix dQ, I would like to find the matrix dU such that dQ = U' * dU + dU' * U. The equation is similar to the ones of this package but does not seem to match any, or does it ?

[baggepinnen]

arec solves A'X + XA - (XB+S)R^(-1)(B'X+S') + Q = 0
so with

X = dU
A = U
Q = -dQ
R = I
B = S = 0

you have a match? Not sure if B=0 is accepted though?

[baggepinnen]

And sylvc solves AX + XB = C, so with A = U', B = U maybe this would work better?

[blegat]

That's what I thought but the issue is that they compute a PSD solution X and it won't work, e.g., if A has eigenvalues of positive real part and Q is PSD (as the system is unstable and a feasible PSD X would give a Lyapunov function).
Here, I am looking for a nonsymmetric X.

[andreasvarga]

This equation has a Lyapunov-like form

A^TX+X^TA = C

and there is no solver for it in MatrixEquations.

However, an approach to solve this equation has been proposed in

H. W. Braden, The Equation $A\sp{T}X-X\sp{T}A=B$, Siam J. Matrix Anal. Appl. 20, 395-402, 1998. [Siam J. Matrix Anal. Appl.]

A generalization of this equation is the Sylvester-like equation

A^TX+X^TB = C

for which a solution method is described in

F. De Teran and Froilan M. Dopico. Consistency and efficient solution of the Sylvester
equation for ⋆-congruence. Electronic J. of Linear Algebra, 22:849–863, 2011.

Interesting information on these equations is also provided in

COMPUTATIONAL METHODS FOR LINEAR MATRIX EQUATIONS
V. SIMONCINI
SIAM ReviewVol. 58, Iss. 3 (2016)10.1137/130912839

[andreasvarga]

For your case, a solution (non-unique) can be determined in two steps. If we denote Y = U^TdU, then we solve for Y first
Y + Y^T = dQ
and then compute the solution by solving
Y = U^TdU
``
In Julia, this can be done in one line

dU = U'\(0.5*Diagonal(diag(dQ))+triu(dQ,1))

[blegat]

Thank you for the detailed answer.
I just thought that if U is upper triangular and we want to constraint dU to be upper triangular as well, there is a unique solution.
It can even be computed in-place by modifying dQ. See jump-dev/DiffOpt.jl#232
I'll go over the reference to see if they mention that upper triangular case (this is in the context of the revision of https://arxiv.org/abs/2206.06135)
Indeed, the one line would work but it would not provide an upper triangular solution.

[andreasvarga]

The general solution according to

Explicit solution of the operator equation A∗X + X∗A=B
Dragan S. Djordjevic
Journal of Computational and Applied Mathematics 200 (2007) 701–704

is

dU = U'\(0.5*dQ)+ Z*U

where Z is any skew-symmetric matrix. I wonder if there exists a skew-symmetric Z to annihilate just the lower triangular part of any given solution.

[andreasvarga]

I verified dU_from_dQ!(dQ, U). It works as expected. Thanks for bringing this to my attention.

[andreasvarga]

I implemented recently two solvers for Lyapunov-like equations. I also implemented a collection of experimental solvers based on Kronecker product expansions. A new patch release will include these solvers.

[blegat]

These solve my problem, thanks!

[andreasvarga]

Given a Cholesky factorization Q = U'U and a symmetric matrix dQ, I would like to find the matrix dU such that dQ = U' * dU + dU' * U.

For your problem, with U upper triangular and nonsigular, a more efficient solution method, which also ensures that dU is upper triangular, is the following:

using LinearAlgebra
n = 3
dQ = Symmetric(rand(n,n))
U = UpperTriangular(rand(n,n))
dU = 0.5*(U'\dQ)   # determine a particular solution
dUl = tril(dU/U,-1)
z = -dUl+dUl'      # determine a skew-symmetric matrix which makes dU+z*U upper triangular
dU1 = triu(dU+z*U)
norm(U'*dU1+dU1'*U-dQ)

[andreasvarga]

Thank you for the detailed answer.
I just thought that if U is upper triangular and we want to constraint dU to be upper triangular as well, there is a unique solution.
It can even be computed in-place by modifying dQ. See jump-dev/DiffOpt.jl#232
I'll go over the reference to see if they mention that upper triangular case (this is in the context of the revision of https://arxiv.org/abs/2206.06135)
Indeed, the one line would work but it would not provide an upper triangular solution.

I wonder if there is a working scheme to determine an upper triangular dU for the case when U is upper triangular but singular.