joint smooth matrix prediction

Joints reconstruct a n by p, q, etc matrices or predictors. The reconstruction can be regularized. # norm( x - u_x v_x^t ) + norm( y - u_y v_y^t) + .... etc

jointSmoothMatrixReconstruction(
  x,
  nvecs,
  parameters,
  iterations = 10,
  subIterations = 5,
  gamma = 1e-06,
  sparsenessQuantile = 0.5,
  positivity = FALSE,
  smoothingMatrix = NA,
  rowWeights = NA,
  repeatedMeasures = NA,
  doOrth = FALSE,
  verbose = FALSE
)

Arguments

x

input list of matrices to be jointly predicted.

nvecs

number of basis vectors to compute

parameters

should be a ncomparisons by 3 matrix where the first two columns define the pair to be matched and the last column defines the weight in the objective function.

iterations

number of gradient descent iterations

subIterations

number of gradient descent iterations in sub-algorithm

gamma

step size for gradient descent

sparsenessQuantile

quantile to control sparseness - higher is sparser

positivity

restrict to positive or negative solution (beta) weights. choices are positive, negative or either as expressed as a string.

smoothingMatrix

a list containing smoothing matrices of the same length as x.

rowWeights

vectors of weights with size n (assumes diagonal covariance)

repeatedMeasures

list of repeated measurement identifiers. this will allow estimates of per identifier intercept.

doOrth

boolean enforce gram-schmidt orthonormality

verbose

boolean option

Value

matrix list each of size p by k is output

Author

Avants BB

Examples

if (FALSE) { # \dontrun{
mat <- replicate(100, rnorm(20))
mat2 <- replicate(100, rnorm(20))
mat <- scale(mat)
mat2 <- scale(mat2)
wt <- 0.666
mat3 <- mat * wt + mat2 * (1 - wt)
params <- matrix(nrow = 2, ncol = 3)
params[1, ] <- c(1, 2, 1)
params[2, ] <- c(2, 1, 1)
x <- list((mat), (mat3))
jj <- jointSmoothMatrixReconstruction(x, 2, params,
  gamma = 1e-4, sparsenessQuantile = 0.5, iterations = 10,
  smoothingMatrix = list(NA, NA), verbose = TRUE
)

# latent feature
mask <- getMask(antsImageRead(getANTsRData("r16")))
spatmat <- t(imageDomainToSpatialMatrix(mask, mask))
smoomat <- knnSmoothingMatrix(spatmat, k = 27, sigma = 120.0)
lfeats <- t(replicate(100, rnorm(3)))
# map these - via matrix - to observed features
n <- sum(mask)
ofeats1 <- (lfeats + rnorm(length(lfeats), 0.0, 1.0)) %*% rbind(rnorm(n), rnorm(n), rnorm(n))
ofeats2 <- (lfeats + rnorm(length(lfeats), 0.0, 1.0)) %*% rbind(rnorm(n), rnorm(n), rnorm(n))
# only half of the matrix contains relevant data
ofeats1[, 1:round(n / 2)] <- matrix(rnorm(round(n / 2) * 100), nrow = 100)
ofeats2[, 1:round(n / 2)] <- matrix(rnorm(round(n / 2) * 100), nrow = 100)
x <- list((ofeats1), (ofeats2))
jj <- jointSmoothMatrixReconstruction(x, 2, params,
  gamma = 0.0001, sparsenessQuantile = 0.75, iterations = 19,
  subIterations = 11,
  smoothingMatrix = list(smoomat, smoomat), verbose = TRUE
)
p1 <- ofeats2 %*% jj$v[[2]]
p2 <- ofeats1 %*% jj$v[[1]]
cor(p1, lfeats)
cor(p2, lfeats)
print(cor(rowMeans(ofeats1), lfeats))
print(cor(rowMeans(ofeats2), lfeats))

# a 2nd example with 3 modalities
imageIDs <- c("r16", "r27", "r30", "r62", "r64", "r85")
images <- list()
feature1Images <- list()
feature2Images <- list()
feature3Images <- list()
ref <- antsImageRead(getANTsRData("r16"))
for (i in 1:length(imageIDs))
{
  cat("Processing image", imageIDs[i], "\n")
  tar <- antsRegistration(ref, antsImageRead(getANTsRData(imageIDs[i])),
    typeofTransform = "Affine"
  )$warpedmov
  images[[i]] <- tar
  feature1Images[[i]] <- iMath(images[[i]], "Grad", 1.0)
  feature2Images[[i]] <- iMath(images[[i]], "Laplacian", 1.0)
  feature3Images[[i]] <- reflectImage(tar, axis = 0, tx = "Affine")$warpedmovout
}
i <- 1
mask <- getMask(antsImageRead(getANTsRData(imageIDs[i])))
mask2 <- iMath(mask, "ME", 2)
spatmat <- t(imageDomainToSpatialMatrix(mask, mask))
smoomat <- knnSmoothingMatrix(spatmat, k = 125, sigma = 100.0)
spatmat2 <- t(imageDomainToSpatialMatrix(mask2, mask2))
smoomat2 <- knnSmoothingMatrix(spatmat2, k = 125, sigma = 100.0)
params <- matrix(nrow = 6, ncol = 3)
params[1, ] <- c(1, 2, 1)
params[2, ] <- c(2, 1, 1)
params[3, ] <- c(1, 3, 1)
params[4, ] <- c(3, 1, 1)
params[5, ] <- c(2, 3, 1)
params[6, ] <- c(3, 2, 1)
mat <- imageListToMatrix(feature1Images, mask)
mat2 <- imageListToMatrix(feature2Images, mask2)
mat3 <- imageListToMatrix(feature3Images, mask)
scl <- F
x <- list(scale(mat, scale = scl), scale(mat2, scale = scl), scale(mat3, scale = scl))
slist <- list(smoomat2, smoomat, smoomat, smoomat, smoomat, smoomat2)

jj <- jointSmoothMatrixReconstruction(x, 4, params,
  positivity = T,
  gamma = 1e-6, sparsenessQuantile = 0.9, iterations = 10,
  smoothingMatrix = slist, verbose = TRUE
)
mm <- makeImage(mask, abs(jj$v[[2]][, 1])) %>% iMath("Normalize")
plot(ref, mm, doCropping = FALSE, window.overlay = c(0.1, 1))

p1 <- mat2 %*% jj$v[[1]]
p2 <- mat %*% jj$v[[2]]
diag(cor(p1, p2))

p1 <- mat3 %*% jj$v[[5]]
p2 <- mat2 %*% jj$v[[6]]
diag(cor(p1, p2))
} # }