Nothing
#' Trace Ratio Optimation
#'
#' This function provides several algorithms to solve the following problem
#' \deqn{\textrm{max} \frac{tr(V^\top A V)}{tr(V^\top B V)} \textrm{ such that } V^\top C V = I}
#' where \eqn{V} is a projection matrix, i.e., \eqn{V^\top V = I}. Trace ratio optimization
#' is pertained to various linear dimension reduction methods. It should be noted that
#' when \eqn{C = I}, the above problem is often reformulated as a generalized eigenvalue problem
#' since it's an easier proxy with faster computation.
#'
#' @param A a \eqn{(p\times p)} symmetric matrix in the numerator term.
#' @param B a \eqn{(p\times p)} symmetric matrix in the denomiator term.
#' @param C a \eqn{(p\times p)} symmetric constraint matrix. If not provided, it is set as identical matrix automatically.
#' @param dim an integer for target dimension. It can be considered as the number of loadings.
#' @param method the name of algorithm to be used. Default is \code{2003Guo}.
#' @param maxiter maximum number of iterations to be performed.
#' @param eps stopping criterion for iterative algorithms.
#'
#' @return a named list containing
#' \describe{
#' \item{V}{a \eqn{(p\times dim)} projection matrix.}
#' \item{tr.val}{an attained maximum scalar value.}
#' }
#'
#' @examples
#' ## simple test
#' # problem setting
#' p = 5
#' mydim = 2
#' A = matrix(rnorm(p^2),nrow=p); A=A%*%t(A)
#' B = matrix(runif(p^2),nrow=p); B=B%*%t(B)
#' C = diag(p)
#'
#' # approximate solution via determinant ratio problem formulation
#' eigAB = eigen(solve(B,A))
#' V = eigAB$vectors[,1:mydim]
#' eigval = sum(diag(t(V)%*%A%*%V))/sum(diag(t(V)%*%B%*%V))
#'
#' # solve using 4 algorithms
#' m12 = trio(A,B,dim=mydim, method="2012Ngo")
#' m09 = trio(A,B,dim=mydim, method="2009Jia")
#' m07 = trio(A,B,dim=mydim, method="2007Wang")
#' m03 = trio(A,B,dim=mydim, method="2003Guo")
#'
#' # print the results
#' line1 = '* Evaluation of the cost function'
#' line2 = paste("* approx. via determinant : ",eigval,sep="")
#' line3 = paste("* trio by 2012Ngo : ",m12$tr.val, sep="")
#' line4 = paste("* trio by 2009Jia : ",m09$tr.val, sep="")
#' line5 = paste("* trio by 2007Wang : ",m07$tr.val, sep="")
#' line6 = paste("* trio by 2003Guo : ",m03$tr.val, sep="")
#' cat(line1,"\n",line2,"\n",line3,"\n",line4,"\n",line5,"\n",line6)
#'
#' @references
#' \insertRef{guo_generalized_2003}{maotai}
#'
#' \insertRef{wang_trace_2007}{maotai}
#'
#' \insertRef{yangqingjia_trace_2009}{maotai}
#'
#' \insertRef{ngo_trace_2012}{maotai}
#'
#' @export
trio <- function(A, B, C, dim=2, method=c("2003Guo","2007Wang","2009Jia","2012Ngo"), maxiter=1000, eps=1e-10){
###################################################################
# not completed yet.
if (missing(C)){
C = diag(nrow(A))
myflag = TRUE
} else {
myflag = FALSE
}
if (!check_symm(A)){
stop("* trio : an input matrix 'A' should be a square, symmetric matrix.")
}
if (!check_symm(B)){
stop("* trio : an input matrix 'B' should be a square, symmetric matrix.")
}
if (!check_symm(C)){
stop("* trio : an input matrix 'C' should be a square, symmetric matrix.")
}
sizes = rep(0,3)
sizes[1] = nrow(A)
sizes[2] = nrow(B)
sizes[3] = nrow(C)
if (length(unique(sizes))!=1){
stop("* trio : all input matrices should be of same size.")
}
if (!myflag){
eigC = eigen(C)
Cinv2 = eigC$vectors%*%diag(1/sqrt(eigC$values))%*%t(eigC$vectors)
A = Cinv2%*%A%*%Cinv2
B = Cinv2%*%B%*%Cinv2
}
# 2009 Jia's note : B should have rank >= (m-d)
if (as.integer(Matrix::rankMatrix(B))<(nrow(B)-dim)){
warning("* trio : null space of 'B' is excessive. trace ratio value may diverge.")
}
###################################################################
# switch case
V = switch(method,
"2007Wang" = trio2007Wang(A, B, dim, eps, maxiter),
"2003Guo" = trio2003Guo(A, B, dim, eps, maxiter),
"2009Jia" = trio2009Jia(A, B, dim, eps, maxiter),
"2012Ngo" = trio2012Ngo(A, B, dim, eps, maxiter))
output = list()
output$V = V
output$tr.val = sum(diag(t(V)%*%A%*%V))/sum(diag(t(V)%*%B%*%V))
return(output)
}
# subroutines -------------------------------------------------------------
#' 2003 Guo et al.
#' Title : A generalized Foley-Sammon transform based on generalized fisher discriminant ...
#'
#' @keywords internal
#' @noRd
trio2003Guo <- function(A, B, dim, eps, maxiter){
## translate into the language
d = dim
Sp = A
Sl = B
## bisection
# 1. initialization
lbd1 = 0; f1 = evalGuoDiff(lbd1, Sp, Sl, d)
lbd2 = 1; f2 = evalGuoDiff(lbd2, Sp, Sl, d)
if (f2 >= 0){
while (f2 > 0){
lbd1 = lbd2; f1 = f2;
lbd2 = lbd2*2; f2 = evalGuoDiff(lbd2, Sp, Sl, d)
}
}
for (i in 1:maxiter){
lbdm = (lbd1+lbd2)/2
fm = evalGuoDiff(lbdm, Sp, Sl, d)
if (fm > 0){
lbd1 = lbdm
f1 = fm
} else {
lbd2 = lbdm
f2 = fm
}
if (abs(lbd1-lbd2) < eps){
break
}
}
lbdm = (lbd1+lbd2)/2
# W = eigen(Sp-lbdm*Sl)$vectors[,1:d] ## use RSpectra for only top 'd' ones
W = RSpectra::eigs(Sp-lbdm*Sl,d,which="LR")$vectors
## let's try to return !
return(W)
}
#' @keywords internal
#' @noRd
evalGuoDiff <- function(lbd, A, B, dim){
W = RSpectra::eigs(A-lbd*B,dim,which="LR")$vectors
# W = eigen(A-lbd*B)$vectors[,1:dim] ## use RSpectra for only top 'd' ones
return(sum(diag(t(W)%*%(A-lbd*B)%*%W)))
}
#' 2007 Wang et al.
#' Title : Trace Ratio vs. Ratio Trace for Dimensionality Reduction
#'
#' @keywords internal
#' @noRd
trio2007Wang <- function(A, B, dim, eps, maxiter){
## translate into this language
Sp = A
St = A+B
m = nrow(A)
d = dim
eigSt = base::eigen(St, symmetric = TRUE)
mm = sum(eigSt$values > 0)
if (mm < 1){
stop("* (A+B) has at least one nonnegative eigenvalues.")
}
U = eigSt$vectors[,1:mm]
## transform into the problem of V now.
Spu = t(U)%*%Sp%*%U
Stu = t(U)%*%St%*%U
Vold = qr.Q(qr(matrix(rnorm(mm*d),ncol=d)))
## main computation
V = cppsub_2007Wang(Vold, mm, d, Spu, Stu, maxiter, eps)
## adjust back to the original problem
W = U%*%V
## let's try to return !
return(W)
}
#' 2009 Jia et al
#' Title : Trace Ratio Problem Revisited (DNM)
#'
#' @keywords internal
#' @noRd
trio2009Jia <- function(A, B, dim, eps, maxiter){
## translate into the language
d = dim
Sp = A
Sl = B
## Decomposed Newton Method
lbdold = 0
for (i in 1:maxiter){
## 1. compute eigendecomposition
eigS = RSpectra::eigs(Sp-lbdold*Sl,d,which="LR")
top.val = eigS$values # top 'd' eigenvalues
top.vec = eigS$vectors # top 'd' eigenvectors
## 2. lbdnew
lbdnew = solve2009Jia(lbdold, top.val, top.vec, Sl)
inc = abs(lbdnew-lbdold)
## 3. updating information
lbdold = lbdnew
if (inc<eps){
break
}
}
## now lbdold is our solution
W = RSpectra::eigs(Sp-lbdold*Sl,d,which="LR")$vectors
## let's try to return !
return(W)
}
#' @keywords internal
#' @noRd
solve2009Jia <- function(lbd, eval, evec, Sl){
d = length(eval)
primesum = 0
for (i in 1:d){
tgt = as.vector(evec[,i])
primesum = primesum - sum(as.vector(Sl%*%tgt)*tgt)
}
sol = (lbd*primesum - sum(eval))/primesum
return(sol)
}
#' 2012 Ngo et al
#' Title : The Trace Ratio Optimization Problem (Newton-Lanczos Algorithm)
#'
#'
#' @keywords internal
#' @noRd
trio2012Ngo <- function(A, B, dim, eps, maxiter){
## in the the language
n = nrow(A)
p = dim
## prepare the initializer
Vold = qr.Q(qr(matrix(rnorm(n*p),ncol=p)))
rhoold = 0
for (i in 1:maxiter){
Vnew = RSpectra::eigs(A-rhoold*B,p,which="LR")$vectors
rhonew = sum(diag(t(Vnew)%*%A%*%Vnew))/sum(diag(t(Vnew)%*%B%*%Vnew))
rhoinc = abs(rhonew-rhoold)
Vold = Vnew
rhoold = rhonew
if (rhoinc < eps){
break
}
}
## let's try to return !
return(Vold)
}
# #
# # simple test
# p = 100
# mydim = 10
# A = matrix(rnorm(p^2),nrow=p); A=A%*%t(A)
# B = matrix(runif(p^2),nrow=p); B=B%*%t(B)
# C = diag(p)
#
#
# library(geigen)
# eigAB = eigen(solve(B,A)) ## geigen(B,A) # we need largest so be careful of the order
# mylist = list()
# V = eigAB$vectors[,1:mydim]; mylist$V = V
# myval = sum(diag(t(V)%*%A%*%V))/sum(diag(t(V)%*%B%*%V)); mylist$tr.val = myval
#
# m12 = trio(A,B,dim=mydim, method="2012Ngo")
# m9 = trio(A,B,dim=mydim, method="2009Jia")
# m7 = trio(A,B,dim=mydim, method="2007Wang")
# m3 = trio(A,B,dim=mydim, method="2003Guo")
#
# # try rstiefel's optimization
# library(rstiefel)
# f1 = function(w){
# trA = sum(diag(t(w)%*%A%*%w)); trB = sum(diag(t(w)%*%B%*%w));
# return(-trA/trB)
# }
# df1 = function(w){
# trA = sum(diag(t(w)%*%A%*%w)); trB = sum(diag(t(w)%*%B%*%w));
# t1 = -2*(A%*%w)*trB + trA*2*(B%*%w);
# t2 = trB^2;
# return(t1/t2)
# }
# f2 = function(w){
# trA = sum(diag(t(w)%*%A%*%w)); trB = sum(diag(t(w)%*%B%*%w));
# return(trB/trA)
# }
# df2 = function(w){
# trA = sum(diag(t(w)%*%A%*%w)); trB = sum(diag(t(w)%*%B%*%w));
# t1 = 2*trA*(B%*%w) - 2*trB*(A%*%w);
# t2 = trA^2;
# return(t1/t2)
# }
# V1 <- optStiefel(f1, df1, Vinit=rustiefel(p,mydim), maxIters = 999, tol=1e-10, verbose=TRUE)
# V2 <- optStiefel(f2, df2, Vinit=rustiefel(p,mydim), maxIters = 999, tol=1e-10, verbose=TRUE)
# print(sprintf("result of (2) inverse : %f", 1/f2(V2)))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.