cc_svd: Conjugate complement of span(X) in span(Z) with respect to...

In gmonette/spida2: Collection of tools developed for the Summer Programme in Data Analysis 2000-2012

Conjugate complement of span(X) in span(Z) with respect to inner product ip

Description

Conjugate complement of span(X) in span(Z) with respect to inner product ip

Usage

cc_svd(X, Z = NULL, ip = NULL, tol = 1e-14)

cc_qr(X, Z = diag(NROW(X)), ip = diag(NROW(X)), tol = 1e-14)

Arguments

X	matrix defining space whose complement is computed. Not necessarily of full column rank
Z	matrix defining space within which the complement is computed. Should be of full column rank. Default: diag(nrow(X))
ip	positive definite matrix defining inner product with respect to which complement is computed. Default: diag(nrow(X))
tol	tolerance (default 1e-14)

Functions

• cc_qr(): Conjugate complement using QR algorithm

Examples

#' cc_svd(cbind(1:3))
cc_svd(cbind(1:3), cbind(1:3,1))
cc_svd(cbind(1:3), cbind(3*c(1,-2,1), 1:3,1))
# challenge with large numbers
K <- cc_svd(cbind(1,1:10,(1:10)^2))
K
ret <- list()
for(i in 1:20){
  cbind(1,(1:10-10^i)^2,1:10) %>% 
    { list(svd = ccor(cc_svd(.), K) -1, qr =ccor(cc_qr(.), K) -1)
    }  %>% lapply(abs) %>% 
    {list(med=lapply(.,median), ave = lapply(.,mean))} %>% unlist->ret[[i]]
}
ret <- do.call(rbind, ret)
ret
library(lattice)
library(latticeExtra)
xyplot(ts(log2(ret)),  scales = list(y = 'same'), ylab = 'cancorr: log base 2 of deviation from 1')
#
# challenge with small numbers
#
i <- 1
K <- cc_svd(cbind(1,1:10,(1:10 == 1)* 10^{-i}))
K
ret <- list()
for(i in 1:30){
  cbind(1,1:10,(1:10 == 1)* 10^{-i}) %>% 
    { list(svd = ccor(cc_svd(.), K) -1, qr =ccor(cc_qr(.), K) -1)
    }  %>% lapply(abs) %>% 
    {list(med=lapply(.,median), ave = lapply(.,mean))} %>% unlist->ret[[i]]
}
ret <- do.call(rbind, ret)
ret

xyplot(ts(log2(ret)),  scales = list(y = 'same'), ylab = 'cancorr: log base 2 of deviation from 1')

#
# challenge with small numbers added to 1
#
i <- 1
K <- cc_svd(cbind(1,1:10,1+(1:10 == 1)* 10^{-i}, (1:10 != 1)))
K
ret <- list()
for(i in 1:30){
  cbind(1,1:10, 1 + (1:10 == 1)* 10^{-i}) %>% 
    { list(svd = ccor(cc_svd(.), K) -1, qr =ccor(cc_qr(.), K) -1)
    }  %>% lapply(abs) %>% 
    {list(med=lapply(.,median), ave = lapply(.,mean))} %>% unlist->ret[[i]]
}
ret <- do.call(rbind, ret)
ret

gmonette/spida2 documentation built on June 12, 2025, 9:44 p.m.