sim_diag: Simultaneously diagonalize two symmetric matrices (1 being...

In aaronjfisher/qp1qc: Solves a possibly non-convex quadratic program with 1 quadratic constraint

Description

Simultaneously diagonalize two symmetric matrices (1 being positive definite)

Usage

sim_diag(M1, M2, eigenM1 = NULL, eigenM2 = NULL, tol = 10^-4,
  return_diags = FALSE)

Arguments

M1	a positive definite symmetric matrix
M2	a symmetric matrix
eigenM1	(optional) the value of eigen(M1), if precomputed and available
eigenM2	(optional) the value of eigen(M2), if precomputed and available
tol	used for error checks
return_diags	should only the diagonalizing matrix be returned.

Details

This function determines an invertible matrix Q such that (Q' M1 Q) is an identity matrix, and (Q' M2 Q) is diagonal, where Q' denotes the transpose of Q. Note, Q is not necessarily orthonormal or symmetric.

Value

If return_diags = FALSE, the matrix Q is returned. Otherwise, a list with the following elements is returned

• diagonalizer - the matrix Q

• inverse_diagonalizer - the inverse of Q

• M1diag - the diagonal elements of (Q' M1 Q)

• M2diag - the diagonal elements of (Q' M2 Q)

Examples

p <- 5
M1 <- diag(p)+2
M2 <- diag(p) + crossprod(matrix(rnorm(p^2),p,p))
M2[1,] <- M2[,1] <- 0
sdb <- sim_diag(M1=M1,M2=M2,tol=10^-10,return_diags=TRUE)
Q <- sdb$diagonalizer
QM1Q <- t(Q) %*% M1 %*% Q
QM2Q <- t(Q) %*% M2 %*% Q
range(QM1Q -diag(sdb$M1diag))
range(QM2Q -diag(sdb$M2diag))
range(Q %*% sdb$inverse_diagonalizer - diag(p))
range(sdb$inverse_diagonalizer %*% Q - diag(p))

# if M1 is not p.d., a warning is produced, but the computation
# proceeds by switching M1 & M2, and then switching them back.
sdb <- sim_diag(M1=M2,M2=M1,tol=10^-10,return_diags=TRUE)

aaronjfisher/qp1qc documentation built on Aug. 22, 2020, 2:35 p.m.