-
Notifications
You must be signed in to change notification settings - Fork 10
/
linear_OLDA.R
executable file
·147 lines (140 loc) · 4.43 KB
/
linear_OLDA.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
#' Orthogonal Linear Discriminant Analysis
#'
#' Orthogonal LDA (OLDA) is an extension of classical LDA where the discriminant vectors are
#' orthogonal to each other.
#'
#' @param X an \eqn{(n\times p)} matrix or data frame whose rows are observations
#' and columns represent independent variables.
#' @param label a length-\eqn{n} vector of data class labels.
#' @param ndim an integer-valued target dimension.
#' @param preprocess an additional option for preprocessing the data.
#' Default is "center". See also \code{\link{aux.preprocess}} for more details.
#'
#' @return a named list containing
#' \describe{
#' \item{Y}{an \eqn{(n\times ndim)} matrix whose rows are embedded observations.}
#' \item{trfinfo}{a list containing information for out-of-sample prediction.}
#' \item{projection}{a \eqn{(p\times ndim)} whose columns are basis for projection.}
#' }
#'
#' @examples
#' ## use iris data
#' data(iris)
#' set.seed(100)
#' subid = sample(1:150, 50)
#' X = as.matrix(iris[subid,1:4])
#' label = as.factor(iris[subid,5])
#'
#' ## compare with LDA
#' out1 = do.lda(X, label)
#' out2 = do.olda(X, label)
#'
#' ## visualize
#' opar <- par(no.readonly=TRUE)
#' par(mfrow=c(1,2))
#' plot(out1$Y, pch=19, col=label, main="LDA")
#' plot(out2$Y, pch=19, col=label, main="Orthogonal LDA")
#' par(opar)
#'
#' @references
#' \insertRef{ye_characterization_2005}{Rdimtools}
#'
#' @rdname linear_OLDA
#' @author Kisung You
#' @concept linear_methods
#' @export
do.olda <- function(X, label, ndim=2, preprocess=c("center","scale","cscale","whiten","decorrelate")){
#------------------------------------------------------------------------
## PREPROCESSING
# 1. data matrix
aux.typecheck(X)
n = nrow(X)
p = ncol(X)
# 2. label vector
label = check_label(label, n)
ulabel = unique(label)
K = length(ulabel)
if (K==1){
stop("* do.olda : 'label' should have at least 2 unique labelings.")
}
if (K==n){
stop("* do.olda : given 'label' has all unique elements.")
}
if (any(is.na(label))||(any(is.infinite(label)))){
stop("* Supervised Learning : any element of 'label' as NA or Inf will simply be considered as a class, not missing entries.")
}
# 3. ndim
ndim = as.integer(ndim)
if (!check_ndim(ndim,p)){
stop("* do.olda : 'ndim' is a positive integer in [1,#(covariates)].")
}
# 4. preprocess
if (missing(preprocess)){
algpreprocess = "center"
} else {
algpreprocess = match.arg(preprocess)
}
#------------------------------------------------------------------------
## COMPUTATION : PRELIMINARY
# 1. preprocessing of data : note that output pX still has (n-by-p) format
tmplist = aux.preprocess.hidden(X,type=algpreprocess,algtype="linear")
trfinfo = tmplist$info
pX = tmplist$pX
# 2. Sb and Sw
# 2-1. Sb
mu_overall = colMeans(pX)
Sb = array(0,c(p,p))
for (i in 1:K){
idxnow = which(label==ulabel[i])
Pi = length(idxnow)/n
mdiff = (colMeans(pX[idxnow,])-mu_overall)
Sb = Sb + Pi*outer(mdiff,mdiff)
}
# 2-2. Sw
Sw = array(0,c(p,p))
for (i in 1:K){
idxnow = which(label==ulabel[i])
Pi = length(idxnow)/n
Si = array(0,c(p,p))
mu_K = colMeans(pX[idxnow,])
for (j in 1:length(idxnow)){
mdiff = (as.vector(pX[idxnow[j],])-mu_K)
Si = Si + outer(mdiff,mdiff)
}
Sw = Sw + Pi*Si
}
# 2-3. pseudo-inverse for Sw; using my function
invSw = aux.pinv(Sw)
#------------------------------------------------------------------------
## COMPUTATION : MAIN PART FOR UNCORRELATED LDA
# 1. initialize
Dmat = as.vector(aux.geigen(Sb, Sw, 1, maximal=TRUE))
# 2. step-by-step computation
diagIp = diag(p)
for (i in 2:ndim){
# 2-1.
if (i==2){
Dmat = matrix(Dmat)
}
tmpA = t(Dmat)%*%invSw%*%Dmat
tmpB = t(Dmat)%*%invSw
# 2-2. solve intermediate inverse problem
tmpsolve = aux.bicgstab(tmpA, tmpB, verbose=FALSE)$x
Pmat = diagIp - (Dmat)%*%tmpsolve
# 2-3. cost function for outer generalized eigenvalue problem and solve
csolution = as.vector(aux.geigen(Pmat%*%Sb, Sw, 1, maximal=TRUE))
Dmat = cbind(Dmat, csolution)
}
# 2-4. remove column names
colnames(Dmat)=NULL
#------------------------------------------------------------------------
## RETURN
# 1. adjust with orthogonalization
projection = aux.adjqr(Dmat)
# 2. return output
result = list()
result$Y = pX%*%projection
result$trfinfo = trfinfo
result$projection = projection
return(result)
}