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 = 0.000001,
  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) {
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 ) )

}