| maxTraceQuotient {lazy.mat} | R Documentation |
tr(t(X) A X) / tr( t(X) B X) t(X) C X = I.Maximizing the Trace Quotient of the form
tr(t(X) A X) / tr( t(X) B X)
subject to t(X) C X = I.
maxTraceQuotient( A, B = diag(nrow(A)), C = diag(nrow(A)), ndim = 1, init = 1, X = NULL, maxiter = 20, eps = 1e-08, epsnu = 1e-08, skipnu = 0, fasteigen = 0, print = 0 ) traceQuotient(X, A, B) quotientTrace(X, A, B) gradtQ(X, A, B) gradQt(X, A, B) Jennrichnu(X, G, C = NULL)
A |
n x n PSD matrix in the numerator |
B |
n x n PSD matrix in the denominator |
C |
n x n PSD matrix to be used as constraints. |
ndim |
# of columns of X matrix |
init |
= 1 to use Generalized Eigenvalue Decomposition |
X |
n x ndim initial value matrix, if any. |
maxiter |
Max # of iterations |
eps |
Convergence criterion for the relative change of the value |
epsnu |
Convergence criterion for nu. |
skipnu |
= 1 to skip calculation of nu in each iteration |
fasteigen |
= 1 to use |
print |
= 1 to print the resuot, = 2 to print the iteration history. |
G |
n x ndim gradient matrix for |
Algorithm 4.1 of Ngo, Bellalij and Saad (2010) is used
with additional C matrix.
Theoretically, when ndim=1, or ndim>1 and B=C,
init=1 will provide the analytic solution using eigen or
generalized eigen decomposition.
Use epsnu<0 to see how the algorithm behaves.
Howerver, the algorithm sometimes does not converge and the value of the
Trace Quotient may increase during the cource of iteration.
This often happens when the C matrix is the same as the B matrix.
Or, the algorithm may oscillate between two values.
The reason is not clear but may be due to negative large eigen values of the G matrix below.
If 0 < epsnu, it is used to check the convergence including the initial value if X.
traceQuotient calculates the Trace Quotient.
quotientTrace calculates the Quotient Trace defined as
tr( inv(t(X) B X) t(X) A X ).
gradtQ calculates the gradient of Trace Quotient
w.r.t X w/o constraints.
gradqT calculates the gradient of Quotient Trace
w.r.t X w/o constraints.
A list of
X Solution matrix
value The value of the trace quotient maximized.
nu Jennrich's nu criterion for convergence.
tXCX The value of the constraint.
conv Convergence code: 1 for convergence,
2 for convergence prior to the iteration,
0 for non-convergence, and -1 for oscillation between two solutions.
Ngo, T. T. , Bellalij, M. , and Saad, Y. (2010) The Trace Ratio Optimization Problem for Dimensionality Reduction. SIAM Journal on Matrix Analysis and Applications. vol. 31. no. 5, pp. 2950-2971.
Jennrich, R. I. (2001) A simple general method for orthogonal rotation. Psychometrika. vol. 66 no. 2, 289-306.
# In the following examples 1, 2, and 3
# q_n_ must be equal to q_n_a, where _n_=1,2,3,
# When ndim=1, q4 must be equal to q4a
# but when ndim > 1, q4 != q4a.
#
# the size of A, B, and C
n=5
set.seed(1701)
A <- matrix(rnorm(n^2), n,n)
A0 <- t(A)%*%A
B <- matrix(rnorm(n^2), n,n)
B0 <- t(B)%*%B
C <- matrix(rnorm(n^2), n,n)
C0 <- t(C)%*%C
rm(A,B,C)
# # of columns of X
ndim <- 1
# problem 1
title <- "q1=x'x / x'x"
A <- diag(n); B <- diag(n)
res1 <- maxTraceQuotient( A, B, ndim=ndim )
x1 <- res1$X; q1 <- res1$value; nu1 <- res1$nu; tXCX1 <- res1$tXCX
cat("\n\n",title, ": ndim =", ndim, "\n")
Print(x1,q1,nu1,tXCX1, fmt="10.6 8.5")
# analytic solution by Eigenvalue Decomposition (GEVD).
res1a <- eigen( A )
x1a <- res1a$vectors[,1:ndim,drop=0]; q1a <- res1a$value[1:ndim]
q1aa <- traceQuotient( x1a, A, B ); q1al <- sum(q1a)/ndim
Print(x1a,q1a, q1aa, q1al, fmt="10.6 8.5")
Print(t(x1a)%*%x1a)
Print("***",round(q1-q1aa,6))
# problem 2
title <- "q2=x'Ax / x'x"
A <- A0; B <- diag(n)
res2 <- maxTraceQuotient( A, B, ndim=ndim )
x2 <- res2$X; q2 <- res2$value; nu2 <- res2$nu; tXCX2 <- res2$tXCX
cat("\n\n",title, ": ndim =", ndim, "\n")
Print(x2,q2,nu2,tXCX2, fmt="10.6 8.5")
# analytic solution by Eigenvalue Decomposition (GEVD).
res2a <- eigen( A )
x2a <- res2a$vectors[,1:ndim]; q2a <- res2a$value[1:ndim]
q2aa <- traceQuotient( x2a, A, B ); q2al <- sum(q2a)/ndim
Print(x2a,q2a, q2aa, q2al, fmt="10.6 8.5")
Print(t(x2a)%*%x2a)
Print("***",round(q2-q2aa,6))
# problem 3
title <- "q3=x'Ax / x'Bx subject to x'Bx=I"
A <- A0; B <- B0
res3 <- maxTraceQuotient( A, B, B, ndim=ndim )
x3 <- res3$X; q3 <- res3$value; nu3 <- res3$nu; tXCX3 <- res3$tXCX
cat("\n\n",title, ": ndim =", ndim, "\n")
Print(x3,q3,nu3,tXCX3, fmt="10.6 8.5")
# analytic solution by Generalized Eigenvalue Decomposition (GEVD).
res3a <- Geigen( A, B, ndim=ndim, check=3 )
x3a <- res3a$vectors; q3a <- res3a$value; tXBX3a <- res3a$tPBP
q3aa <- traceQuotient( x3a, A, B ); q3al <- sum(res3a$values)/ndim
Print(x3a,q3a, q3aa, q3al, fmt="10.6 8.5")
Print(tXBX3a)
Print("***",round(q3-q3aa,6))
# problem 4: Can be solved by Generalized Eigen Problem when ndim=1.
title <- "q4=tr(X'AX) / tr(XB'X) subject to X'CX=I"
A <- A0; B <- B0; C <- C0
res4 <- maxTraceQuotient( A, B, C, ndim=ndim )
x4 <- res4$X; q4 <- res4$value; nu4 <- res4$nu; tXCX4 <- res4$tXCX
cat("\n\n",title, ": ndim =", ndim, "\n")
Print(x4,q4,nu4,tXCX4, fmt="10.6 8.5")
# analytic solution by Eigenvalue Decomposition (GEVD).
res4a <- Geigen( A, B, ndim=ndim, check=4 )
x4a <- res4a$vectors; q4a <- res4a$value; tXBX4a <- res4a$tPBP
# rescale x
x4an <- x4a/c(sqrt(t(x4a)%*%C%*%x4a))
q4aa <- traceQuotient( x4a, A, B ); q4al <- sum(res4a$values)/ndim
Grad <- gradtQ( x4a, A, B )
nuB <- Jennrichnu( x4a, Grad, B )
nuC <- Jennrichnu( x4a, Grad, C )
Print(x4a,x4an)
Print(q4a, q4aa, q4al ,nuB, nuC)
Print(tXBX4a,t(x4a)%*%C%*%x4a, fmt="8.5")
Print("***",round(q4-q4aa,6))
ndim=2
# problem 1
title <- "q1=x'x / x'x"
A <- diag(n); B <- diag(n)
res1 <- maxTraceQuotient( A, B, ndim=ndim )
x1 <- res1$X; q1 <- res1$value; nu1 <- res1$nu; tXCX1 <- res1$tXCX
cat("\n\n",title, ": ndim =", ndim, "\n")
Print(x1,q1,nu1,tXCX1, fmt="10.6 8.5")
# analytic solution by Eigenvalue Decomposition (GEVD).
res1a <- eigen( A )
x1a <- res1a$vectors[,1:ndim,drop=0]; q1a <- res1a$value[1:ndim]
q1aa <- traceQuotient( x1a, A, B ); q1al <- sum(q1a)/ndim
Print(x1a,q1a, q1aa, q1al, fmt="10.6 8.5")
Print(t(x1a)%*%x1a)
Print("***",round(q1-q1aa,6))
# problem 2
title <- "q2=x'Ax / x'x"
A <- A0; B <- diag(n)
res2 <- maxTraceQuotient( A, B, ndim=ndim )
x2 <- res2$X; q2 <- res2$value; nu2 <- res2$nu; tXCX2 <- res2$tXCX
cat("\n\n",title, ": ndim =", ndim, "\n")
Print(x2,q2,nu2,tXCX2, fmt="10.6 8.5")
# analytic solution by Eigenvalue Decomposition (GEVD).
res2a <- eigen( A )
x2a <- res2a$vectors[,1:ndim]; q2a <- res2a$value[1:ndim]
q2aa <- traceQuotient( x2a, A, B ); q2al <- sum(q2a)/ndim
Print(x2a,q2a, q2aa, q2al, fmt="10.6 8.5")
Print(t(x2a)%*%x2a)
Print("***",round(q2-q2aa,6))
# problem 3
title <- "q3=x'Ax / x'Bx subject to x'Bx=I"
A <- A0; B <- B0
res3 <- maxTraceQuotient( A, B, B, ndim=ndim )
x3 <- res3$X; q3 <- res3$value; nu3 <- res3$nu; tXCX3 <- res3$tXCX
cat("\n\n",title, ": ndim =", ndim, "\n")
Print(x3,q3,nu3,tXCX3, fmt="10.6 8.5")
# analytic solution by Generalized Eigenvalue Decomposition (GEVD).
res3a <- Geigen( A, B, ndim=ndim, check=3 )
x3a <- res3a$vectors; q3a <- res3a$value; tXBX3a <- res3a$tPBP
q3aa <- traceQuotient( x3a, A, B ); q3al <- sum(res3a$values)/ndim
Print(x3a,q3a, q3aa, q3al, fmt="10.6 8.5")
Print(tXBX3a)
Print("***",round(q3-q3aa,6))
# problem 4: Can NOT be solved by Generalized Eigen Problem when ndim > 1.
title <- "q4=tr(X'AX) / tr(XB'X) subject to X'CX=I"
A <- A0; B <- B0; C <- C0
res4 <- maxTraceQuotient( A, B, C, ndim=ndim )
x4 <- res4$X; q4 <- res4$value; nu4 <- res4$nu; tXCX4 <- res4$tXCX
cat("\n\n",title, ": ndim =", ndim, "\n")
Print(x4,q4,nu4,tXCX4, fmt="10.6 8.5")
# analytic solution by Eigenvalue Decomposition (GEVD).
res4a <- Geigen( A, B, ndim=ndim, check=4 )
x4a <- res4a$vectors; q4a <- res4a$value; tXBX4a <- res4a$tPBP
q4aa <- traceQuotient( x4a, A, B ); q4al <- sum(res4a$values)/ndim
Grad <- gradtQ( x4a, A, B )
nuB <- Jennrichnu( x4a, Grad, B )
nuC <- Jennrichnu( x4a, Grad, C )
Print(x4a, q4a, q4aa, q4al ,nuB, nuC)
Print(tXBX4a,t(x4a)%*%C%*%x4a, fmt="8.5")
Print("***",round(q4-q4aa,6))