## Raw procedure for learning resonant basis (adaptive filters / generative models) by computing elastic equilibria # # Run in R statistical computing environment with 'source("csca.r")' in proper directory. Use getwd() and setwd(). # With intelligent coupling Q = E{xx^H}^-1 between systems x_i, and with zero response/consuming constants R, # the models phi_i would converge to principal components of data, or at least to the subspace spanned by them. # The algorithm is extended to approximate completely distributed models with simpler diagonal coupling q_i = constant # by consuming the data for activity, R > 0, which separates the models. The data can be complex valued, # and the default constraint Re(x_i) > 0, Im(x_i) = 0 elicits sparse components. Other constraints are possible, # and complex coupling constants would simulate inductive and conductive phenomena in frequency space. # The models summarize steady states of rapid gradient descents of cybernetic objective function # J(x,u) = 1/2 x^H A x - x^H B u, i.e. dx/dt = - G A x + G B u, i.e. xbar = A^-1 B u = phi^H ubar, # where A = (Q^-1 + R E{xbarxbar^H}) and B = E{xbarubar^H}, # and data is consumed for activity, ubar = u - phi R xbar. # Heikki Hyötyniemi & Petri Lievonen -- http://neocybernetics.com/ #library('Matrix') # use Matrix library instead of standard 'matrix' command # for larger sparse data and matrices of specific form # Set: n <- 50 # number of models/features/systems/attractors x_i lambda <- 0.05 # learning rate, near zero: inverse of adaptation time scale in sample time; reduce by several magnitudes if problems complexModels <- TRUE # operate on complex field. stores real and imaginary parts separately instead of using built-in complexSystem <- FALSE # complex data types to illustrate explicit calculations and workaround unimplemented methods. loadData <- TRUE # load data from file centerSamples <- FALSE # center samples to zero mean normalizeSamples <- FALSE # normalize sample vectors (and thus their component sum of squares) to unity whitenSamples <- FALSE # whiten data in its main principal subspace centerVariables <- TRUE # center data dimensions (recommended for numerical stability) normalizeVariables <- FALSE # normalize data dimension variances to unity (recommended for unknown data) initModels <- TRUE # initialize model matrices, set false to continue iteration offlineProcessing <- TRUE # adapt to all data simultaneously; alternatively process only one sample at a time breakSymmetry <- TRUE # enhance sparsity by truncating negative Re(x_i) to zero (rectifying), so that linear subspaces become convex rays searchCoherence <- TRUE # simplify sparsity by constraining Im(x_i) to zero (this constant is not implemented yet) maskA <- FALSE # elicit actual principal components by structural constraints q_i <- 4 # coupling constants' initial values intelligentCoupling = FALSE # Q = E{xx^H}^-1, set response constants r_i to zero to prevent exhaustion regulatedCoupling = TRUE # q_i = 1/E{|x_i|^2} will keep the models near unit length and maximize perceptions E{|x_i|^2} r_i <- 1 # response constants between 0 (decoupled) and 1 (consuming data for activity) (start small if necessary) showModels = TRUE # show models as images showStructure = TRUE # show higher-order structure as images, matrices phi^H phi and E{xx^H} for example showCriterion = TRUE # plot some objective function of iteration, showing convergence for example showEstimate = TRUE # show some data estimates by current models showBalance = TRUE # show some data balances (natural residuals ubar and model activities xbar) by current models showLabels = TRUE # show labels for different plots if (loadData) { m <- 256 # number of data dimensions u_j k <- 500 # number of data samples u # example data from http://neocybernetics.com/research/digits500.txt U_re <- matrix(scan("digits500.txt"), m, k) # load real-valued data samples to columns of matrix Re(U) U_im <- 0 # try library 'matlab' or 'R.matlab' for .mat data files imageW <- imageH <- 16 # image width and height for this special data # also all kinds of transformations are possible, such as generating higher combinatorics, ## U <- mvfft(U_re); U_re <- Re(U); U_im <- Im(U) # or using Fourier-transform, but must be inverse-transformed later #m <- 1024; k <- 8940; imageW <- 32; imageH <- 32 #U_re <- matrix(scan("digits8940.txt"), m+1, k) # example data from http://neocybernetics.com/research/digits8940.zip #U_re <- U_re[1:m,] # last dimension contains labeling data, removed here #U_re <- (U_re+1)/2 # some example data sets courtesy of http://www.cs.nyu.edu/~roweis/data.html #m <- 256; k <- 1100*10; imageW <- 16; imageH <- 16 #library('R.matlab'); mat <- readMat("usps_all.mat"); U_re <- matrix(mat$data,m) / 256 # class labels available, see dim(mat$data) #m <- 1024; k <- 10000 #50000; #imageW <- 32; imageH <- 32 #library('R.matlab'); mat <- readMat("patches.mat"); U_re <- matrix(mat$X[(m*20000+1):(m*20000+m*k)],m,k) #m <- 320; k <- 1404; imageW <- 20; imageH <- 16 #library('R.matlab'); mat <- readMat("binaryalphadigs.mat"); U_re <- matrix(unlist(mat$dat),nrow=m) # extra dimensions available, see summary(mat) #m <- 560; k <- 1965; imageW <- 20; imageH <- 28 #library('R.matlab'); mat <- readMat("frey_rawface.mat"); U_re <- mat$ff #/255 #log(mat$ff) #U_re <- U_re/rowMeans(U_re) - 1; # == (U_re-rowMeans(U_re))/rowMeans(U_re) # cybernetic preprocessing #1 #U_re <- log(U_re); #U_re <- U_re - rowMeans(U_re); # == log(U_re/exp(rowMeans(log(U_re)))) # cybernetic preprocessing #2 #centerSamples <- TRUE; #centerSamples <- TRUE; normalizeSamples <- TRUE; centerVariables <- TRUE; normalizeVariables <- TRUE #m <- 4096; k <- 400; imageW <- 64; imageH <- 64 #library('R.matlab'); mat <- readMat("olivettifaces.mat"); U_re <- mat$faces / 256 #centerSamples <- TRUE; normalizeSamples <- TRUE; centerVariables <- TRUE; normalizeVariables <- TRUE # extra dimensions available, see summary(mat) #m <- 784; k <- 60000; imageW <- 28; imageH <- 28 #library('R.matlab'); mat <- readMat("mnist_all_mod.mat") #U_re <- t( rbind(mat$train0,mat$train1,mat$train2,mat$train3,mat$train4,mat$train5,mat$train6,mat$train7,mat$train8,mat$train9) ) # use sparse Matrix library instead to save memory! #m <- 784; k <- 6006; imageW <- 28; imageH <- 28 #library('R.matlab'); mat <- readMat("mnist_bw_small.mat") #U_re <- t(mat$X.train) #U_re <- U_re[1:m,] # last dimension contains labeling data, not used yet #m <- 10304; k <- 575; imageW <- 112; imageH <- 92 #library('R.matlab'); mat <- readMat("umist_cropped.mat"); U_re <- matrix(mat$facedat[[1]],m) / 256 #for (i in 2:20) { # U_re <- cbind( U_re, matrix(mat$facedat[[i]],m)/256 ) #} #centerSamples <- TRUE; normalizeSamples <- TRUE; centerVariables <- TRUE; normalizeVariables <- TRUE # extra dimensions available, see summary(mat) #m <- 1024; k <- 1428; imageW <- 32; imageH <- 32 #library('R.matlab'); mat <- readMat("lights=27.mat"); U_re <- t(mat$fea) / 256 #centerSamples <- TRUE; normalizeSamples <- TRUE; centerVariables <- TRUE; normalizeVariables <- TRUE # extra dimensions available, see summary(mat) #m <- 1024; k <- 400; imageW <- 32; imageH <- 32 #library('R.matlab'); mat <- readMat("ORL_32x32.mat"); U_re <- t(mat$fea) / 256 #centerSamples <- TRUE; normalizeSamples <- TRUE; centerVariables <- TRUE; normalizeVariables <- TRUE # extra dimensions available, see summary(mat) #m <- 4096; k <- 400; imageW <- 64; imageH <- 64 #library('R.matlab'); mat <- readMat("ORL_64x64.mat"); U_re <- t(mat$fea) / 256 #centerSamples <- TRUE; normalizeSamples <- TRUE; centerVariables <- TRUE; normalizeVariables <- TRUE # extra dimensions available, see summary(mat) #m <- 1024; k <- 11554; imageW <- 32; imageH <- 32 #library('R.matlab'); mat <- readMat("PIE_32x32.mat"); U_re <- t(mat$fea) / 256 #centerSamples <- TRUE; normalizeSamples <- TRUE; centerVariables <- TRUE; normalizeVariables <- TRUE # extra dimensions available, see summary(mat) #m <- 100; k <- 16242; imageW <- 10; imageH <- 10 #library('R.matlab'); mat <- readMat("20news_w100.mat"); U_re <- as.matrix(mat$documents) # data should be preprocessed properly # extra dimensions available, see summary(mat) ##library('R.matlab'); mat <- readMat("grolier15276.mat"); ##library('R.matlab'); mat <- readMat("pnas_all.mat"); ##library('R.matlab'); mat <- readMat("nips12raw_str602.mat"); ##library('R.matlab'); mat <- readMat("nips_1-17.mat"); ##library('R.matlab'); mat <- readMat("visterms_1000.mat"); if (centerSamples) { # optionally center sample vectors (remove the DC component) sampleMeans_re <- colMeans(U_re) # calculate sample means U_re <- t(t(U_re)-sampleMeans_re) # center samples by substracting the means if (complexModels) {} # center also imaginary component if present } if (normalizeSamples) { # optionally normalize sample vectors to unit length sampleSqMeans <- colMeans(U_re^2) # calculate means of sample squares (std.dev^2, RMS^2) sampleSqMeans[sampleSqMeans==0] <- 1 # prevent division by zero U_re <- t(t(U_re)/sqrt(sampleSqMeans)) # normalize sample vectors to unit length if (complexModels) {} # normalize properly with sqrt(Mod()) if imaginary component is present } if (whitenSamples) { # experiment with prewhitening pca <- eigen(tcrossprod(U_re)/k) nbasis <- floor(0.25*length(pca$values)) invvar <- 1/sqrt(pca$values[1:nbasis]) #invvar <- rep(1,nbasis) basis <- pca$vectors[,1:nbasis] #U_re <- basis %*% diag(invvar) %*% t(basis) %*% U_re U_re <- diag(invvar) %*% t(basis) %*% U_re m <- dim(U_re)[1] } maxU <- max(abs(range(U_re))) # store maximum value for drawing if (centerVariables) { # for numerical stability it is better to center the data meanU_re <- rowMeans(U_re) # calculate data dimension averages to set affine subspaces #meanU_re <- exp(rowMeans(log(U_re))) # alternatively, use geometric means for positive data U_re <- U_re - meanU_re # center the data by substracting the mean (store for reversal) if (complexModels) { # center also imaginary component if present #meanU_im <- rowMeans(U_im) #U_im <- U_im - meanU_im } } else { meanU_re <- 0; meanU_im <- 0 } #U_re[m,] <- 1 # experiment with setting some dimension to constant unity to facilitate affine modeling if (normalizeVariables) { # if nothing more specific is known about the data distribution if (complexModels) { varU <- rowMeans(U_re^2 + U_im^2) # calculate data dimension variances } else { varU <- rowMeans(U_re^2) # calculate data dimension variances } #varU <- mean(varU) # optionally treat several dimensions as similar #varU <- meanU_re^2 # cybernetic standardization: divide deviations from means by those means varU[varU==0] <- 1 # prevent division by zero U_re <- U_re/sqrt(varU) # normalize data by dividing by standard deviations if (complexModels) { # normalize properly with sqrt(Mod()) if imaginary component is present U_im <- U_im/sqrt(varU) } } else { varU <- 1 } maxUbar <- max(abs(range(U_re))) # store maximum value for drawing if (complexModels) { # a special example for this data: U_im <- matrix(0,m,k) # set the first 10 imaginary components of data to imply digit U_im[1:10,] <- (row(U_im)%%10 == col(U_im)%%10)[1:10,] meanU_im <- 0 #U_im[,1:(k-1)] <- U_re[,2:k] # second special example: set imaginary component to denote the next real sample #U_im[,1:(k-1)] <- U_re[,2:k] - U_re[,1:(k-1)] # or set imaginary component to denote the gradient ##meanU_im <- meanU_re ##limits <- range(U_re) # third special example: ###U_im <- sin( (U_re - limits[1])/(limits[2]-limits[1]) * pi - pi/2 ) # one can also experiment coding values to ##U_im <- sin( (U_re - limits[1])/(limits[2]-limits[1]) * pi - pi/2 ) * ( (U_re - limits[1])/(limits[2]-limits[1]) * 2 - 1 ) ###U_re <- cos( (U_re - limits[1])/(limits[2]-limits[1]) * pi - pi/2 ) # phases/arguments and confidences to magnitudes/modulus ##U_re <- cos( (U_re - limits[1])/(limits[2]-limits[1]) * pi - pi/2 ) * ( (U_re - limits[1])/(limits[2]-limits[1]) * 2 - 1 ) } else { U_im <- 0 } #U_re[1:10,] <- (row(U_re)%%10 == col(U_re)%%10)[1:10,] # also real dimensions could be used for extra data (this overwrites for simplicity) } # for natural positive data consider instead logarithms and geometric means of data, # for example dividing deviations from means by those means. (see "freyface.mat" above for example) # for time series study vectorizing successive samples together if(initModels) { ExuH_re <- 1/n*matrix(runif(n*m)-0.5,n,m) # initialize matrix Re(E{xu^H}) with small random real numbers (1/sqrt(m), for example) ##ExuH_re <- ExuH_re * t(U_re[,sample.int(k,n)]) # could try different initializations if (complexModels) { ExuH_im <- 1/n*matrix(runif(n*m)-0.5,n,m) # initialize matrix Im(E{xu^H}) with small random imaginary numbers (1/sqrt(m), for example) } else { ExuH_im <- 0 } ExxH_re <- 1/n*matrix(runif(n*n),n,n) # initialize matrix E{xx^H} with small random numbers (or just a diagonal) #invExxH_re <- solve(ExxH_re) if (complexSystem) { ExxH_im <- 1/n*matrix(runif(n*n),n,n) # initialize matrix E{xx^H} with small random numbers (or just a diagonal) } else { ExxH_im <- 0 } Q_re <- q_i*diag(n) # initialize matrix Q with equal coupling constants in diagonal, just for this case if (complexSystem) { Q_im <- 0*diag(n) } else { Q_im <- 0 } if (intelligentCoupling) { # if intelligent coupling is assumed Q_re <- solve(ExxH_re) # initialize coupling matrix as inverse of activity covariance matrix } if (complexSystem) { phi_re <- t(Q_re %*% ExuH_re - Q_im %*% ExuH_im) # initialize real parts of models/filters as columns of matrix Re(phi) phi_im <- -t(Q_re %*% ExuH_im + Q_im %*% ExuH_re) # imaginary part is negative due to complex conjugation ^H } else if (complexModels) { phi_re <- t(Q_re %*% ExuH_re) # initialize real parts of models/filters as columns of matrix Re(phi) phi_im <- -t(Q_re %*% ExuH_im) # initialize imaginary parts of models/filters as columns of matrix Im(phi) } else { phi_re <- t(Q_re %*% ExuH_re) # initialize real parts of models/filters as columns of matrix Re(phi) phi_im <- 0 } R_re <- r_i # initialize matrix R with equal response constants in diagonal (or replace with a scalar) if (complexSystem) { R_im <- 0 } else { R_im <- 0 } Xbar_re <- matrix(0,n,k) # initialize perceptions if (complexSystem) { Xbar_im <- matrix(0,n,k) } else { Xbar_im <- 0 } Ubar_re <- matrix(0,m,k) # initialize effective variables if (complexModels) { Ubar_im <- matrix(0,m,k) } else { Ubar_im <- 0 } } # should make more effort to generate different scenarios and bootstraps, # to reduce the initial transient, for example #jet.colors <- colorRampPalette(c("blue","white","red"), space="Lab") # create custom color palette jet.colors <- colorRampPalette(c("black","white"), space="Lab") # alternative custom color palette bwr <- jet.colors(255) # palette ranges from blue (-) to white (0) to red (+) J <- c(0) # initialize vector for storing some cost criterion to be plotted #Jprev <- sum(diag(ExxH_re)) # Main loop: for (p in 1:10000) { # iterate for some time or when Esc or Ctrl+C pressed currentSample = (p-1)%%k + 1 # index of current sample: iteration number p limited to number of samples k # create/clear a graphics window with 5x5=25 tiled image blocks gridSize = ceiling(sqrt(n*showModels + 2*showStructure + showCriterion + 2*showEstimate + 2*showBalance)) if (complexModels) { par(mfrow=c(gridSize,gridSize*2), mar=c(1,0.1,0.1,0.1)) #edit to fit } else { par(mfrow=c(gridSize,gridSize), mar=c(1,0.1,0.1,0.1)) #edit to fit #par(mfrow=c(6,6), mar=c(1,1,1,1)) #edit to fit } # Filter data u to perceptions xbar: if (offlineProcessing) { if (complexSystem) { # let's calculate the slow way for simplicity varphiH <- Conj(solve(Conj(solve(matrix(complex(real=Q_re,imag=Q_im),n,n))) + complex(real=R_re,imag=R_im) * matrix(complex(real=ExxH_re,imag=ExxH_im),n,n))) %*% matrix(complex(real=ExuH_re,imag=ExuH_im),n,m) Xbar_re <- Re(varphiH) %*% U_re - Im(varphiH) %*% U_im # calculate steady state for each u Xbar_im <- Im(varphiH) %*% U_re + Re(varphiH) %*% U_im # calculate steady state for each u } else if (complexModels) { Xbar_re <- solve(solve(Q_re) + R_re * ExxH_re) %*% (ExuH_re %*% U_re - ExuH_im %*% U_im) # calculate steady state for each u #Xbar_im <- 0 # not calculated: only coherent solutions are searched for #Xbar_im <- solve(solve(Q_re) + R_re * ExxH_re) %*% (ExuH_im %*% U_re + ExuH_re %*% U_im) # calculate steady state for each u #Xbar_re <- (Xbar_re) + Xbar_im # results in rotation, if the next sample is coded in imaginary component #Xbar_re[1:5,] <- Xbar_im[1:5,] # or manipulate only some of the filters } else { Xbar_re <- solve(solve(Q_re) + R_re * ExxH_re) %*% ExuH_re %*% U_re # calculate steady state for each u ##Xbar_re <- solve(diag(n) + R_re %*% phiTphi_re) %*% t(phi_re) %*% U_re # equivalent alternative formulation } } else { # process one sample at a time if (complexSystem) { # let's calculate the slow way for simplicity varphiH <- Conj(solve(Conj(solve(matrix(complex(real=Q_re,imag=Q_im),n,n))) + complex(real=R_re,imag=R_im) * matrix(complex(real=ExxH_re,imag=ExxH_im),n,n))) %*% matrix(complex(real=ExuH_re,imag=ExuH_im),n,m) Xbar_re[,currentSample] <- Re(varphiH) %*% U_re[,currentSample] - Im(varphiH) %*% U_im[,currentSample] # calculate steady state for this u Xbar_im[,currentSample] <- Im(varphiH) %*% U_re[,currentSample] + Re(varphiH) %*% U_im[,currentSample] # calculate steady state for this u } else if (complexModels) { Xbar_re[,currentSample] <- solve(solve(Q_re) + R_re * ExxH_re) %*% (ExuH_re %*% U_re[,currentSample] - ExuH_im %*% U_im[,currentSample]) # calculate steady state for this u } else { Xbar_re[,currentSample] <- solve(solve(Q_re) + R_re * ExxH_re) %*% ExuH_re %*% U_re[,currentSample] # calculate steady state for this u # experiment optimizations by matrix inversion lemma #invAold <- solve(solve(Q_re) + R_re * ExxH_re) #invAnew <- invAold - invAold #Xbar_re[,currentSample] <- invA %*% ExuH_re %*% U_re[,currentSample] # calculate steady state for this u } } # Apply nonlinearity to perceptions xbar: if (breakSymmetry) { # one way to elicit sparse solutions: if (offlineProcessing) { if (complexSystem) { #Xbar_re <- Xbar_re - sqrt(diag(ExxH_re)) #Xbar_re <- Xbar_re - 1/sqrt(diag(Q_re)) #Xbar_re[Xbar_re<(1/sqrt(diag(Q_re)))] <- 0 # truncate negative (or smaller than some constant) inner products to zero, or multiply by small number #Xbar_im[Xbar_re<0] <- 0 # truncate negative (or smaller than some constant) inner products to zero, or multiply by small number Xbar_re[Xbar_re<0] <- 0 # proper nonlinearity should be applied recursively with active models #Xbar_re <- t(phi) %*% U - R_re * phiTphi %*% Xbar_re # but not possible in a linear matrix form #Xbar_re <- matrix(0,n,k) Xbar_im <- matrix(0,n,k) } else if (complexModels) { Xbar_re[Xbar_re<0] <- 0 } else { Xbar_re[Xbar_re<0] <- 0 # one could study applying approximations a few times recursively #Xbar_re <- Xbar_re^2 #Xbar_re <- log(cosh(Xbar_re)) #Xbar_re <- Xbar_re - tanh(Xbar_re) #Xbar_re <- Xbar_re*tanh(Xbar_re) } } else { if (complexSystem) { Xbar_re[Xbar_re[,currentSample]<0,currentSample] <- 0 # truncate negative (or smaller than some constant) inner products to zero Xbar_im[,currentSample] <- 0 } else if (complexModels) { Xbar_re[Xbar_re[,currentSample]<0,currentSample] <- 0 # truncate negative (or smaller than some constant) inner products to zero } else { Xbar_re[Xbar_re[,currentSample]<0,currentSample] <- 0 # truncate negative (or smaller than some constant) inner products to zero } } } # testing: #Xbar_re <- (Xbar_re + Xbar_re*diag(crossprod(phi_re))) * diag(Q_re) #Xbar_re <- (Xbar_re + Xbar_re*diag(crossprod(phi_re))) * 1 #Xbar_re[abs(Xbar_re)<0.3] <- 0 # possible additional sparsifying constraints, based on number of nonzeros, lower quantiles etc. #Xbar_re[abs(Xbar_re)<0.1] <- 0 # possible additional sparsifying constraints, based on number of nonzeros, lower quantiles etc. #Xbar_re[Xbar_re>0.1] <- (Xbar_re/diag(Q_re)/diag(crossprod(phi_re)) + 0.1/diag(Q_re)/diag(crossprod(phi_re)))[Xbar_re>0.1] #Xbar_re[Xbar_re>0.3] <- (Xbar_re/diag(crossprod(phi_re)) + 0.3/diag(crossprod(phi_re)))[Xbar_re>0.3] #Xbar_re[Xbar_re<0.1] <- (Xbar_re/diag(Q_re)/diag(crossprod(phi_re)) - 0.1/diag(Q_re)/diag(crossprod(phi_re)))[Xbar_re<0.1] #Xbar_re[Xbar_re<0.3] <- (Xbar_re/diag(crossprod(phi_re)) - 0.3/diag(crossprod(phi_re)))[Xbar_re<0.3] #Xbar_re[Xbar_re<0.1] <- Xbar_re[Xbar_re<0.1] + (Xbar_re[Xbar_re<0.1]-0.1)/diag(crossprod(phi_re)) #diag(R_re) <- diag(crossprod(phi_re)) #Xbar_re <- Xbar_re^2 # Xbar_re <- sqrt(Xbar_re) # Substract implicit response uest from u to get effective variables ubar: if (offlineProcessing) { if (complexSystem) { if (R_im == 0) { Ubar_re <- U_re - R_re * (phi_re %*% Xbar_re - phi_im %*% Xbar_im) # calculate the effective data by substracting the effect of Ubar_im <- U_im - R_re * (phi_im %*% Xbar_re + phi_re %*% Xbar_im) # steady states: the estimate that is implicit in consuming/dampening response } else { Ubar_re <- U_re - (R_re * (phi_re %*% Xbar_re - phi_im %*% Xbar_im) - R_im * (phi_im %*% Xbar_re + phi_re %*% Xbar_im)) Ubar_im <- U_im - (R_re * (phi_im %*% Xbar_re + phi_re %*% Xbar_im) + R_im * (phi_re %*% Xbar_re - phi_im %*% Xbar_im)) } } else if (complexModels) { Ubar_re <- U_re - R_re * phi_re %*% Xbar_re Ubar_im <- U_im - R_re * phi_im %*% Xbar_re } else { Ubar_re <- U_re - R_re * phi_re %*% Xbar_re } } else { if (complexSystem) { Ubar_re[,currentSample] <- U_re[,currentSample] - R_re * (phi_re %*% Xbar_re[,currentSample] - phi_im %*% Xbar_im[,currentSample]) Ubar_im[,currentSample] <- U_im[,currentSample] - R_re * (phi_im %*% Xbar_re[,currentSample] + phi_re %*% Xbar_im[,currentSample]) } else if (complexModels) { Ubar_re[,currentSample] <- U_re[,currentSample] - R_re * phi_re %*% Xbar_re[,currentSample] Ubar_im[,currentSample] <- U_im[,currentSample] - R_re * phi_im %*% Xbar_re[,currentSample] } else { Ubar_re[,currentSample] <- U_re[,currentSample] - R_re * phi_re %*% Xbar_re[,currentSample] } } # Update system/model matrices essential for filtering and generative modeling: if (offlineProcessing) { # E{xu^H}: if (complexSystem) { ExuH_re <- (1-lambda) * ExuH_re + lambda * (tcrossprod(Xbar_re,Ubar_re) + tcrossprod(Xbar_im,Ubar_im))/k # E{xu^T}_re ~ (Xbar_re %*% t(Ubar_re) - Xbar_im %*% t(-Ubar_im)) / k ExuH_im <- (1-lambda) * ExuH_im + lambda * (tcrossprod(Xbar_im,Ubar_re) - tcrossprod(Xbar_re,Ubar_im))/k # E{xu^T}_im ~ (Xbar_re %*% t(-Ubar_im) + Xbar_im %*% t(Ubar_re)) / k } else if (complexModels) { ExuH_re <- (1-lambda) * ExuH_re + lambda * tcrossprod(Xbar_re,Ubar_re)/k # E{xu^T}_re ~ Xbar_re %*% t(Ubar_re) / k ExuH_im <- (1-lambda) * ExuH_im - lambda * tcrossprod(Xbar_re,Ubar_im)/k # E{xu^T}_im ~ Xbar_re %*% t(-Ubar_im) / k } else { ExuH_re <- (1-lambda) * ExuH_re + lambda * tcrossprod(Xbar_re,Ubar_re)/k # E{xu^T}_re ~ Xbar_re %*% t(Ubar_re) / k } # E{xx^H}: if (complexSystem) { ExxH_re <- (1-lambda) * ExxH_re + lambda * (tcrossprod(Xbar_re) + tcrossprod(Xbar_im))/k # E{xx^T}_re ~ (Xbar_re %*% t(Xbar_re) - Xbar_im %*% t(-Xbar_im)) / k ExxH_im <- (1-lambda) * ExxH_im + lambda * (tcrossprod(Xbar_im,Xbar_re) - tcrossprod(Xbar_re,Xbar_im))/k # E{xx^T}_im ~ (Xbar_re %*% t(-Xbar_im) + Xbar_im %*% t(Xbar_re)) / k } else { ExxH_re <- (1-lambda) * ExxH_re + lambda * tcrossprod(Xbar_re)/k # E{xx^T} ~ Xbar_re %*% t(Xbar_re) / k } } else { # E{xu^H}: if (complexSystem) { ExuH_re <- (1-lambda) * ExuH_re + lambda * (tcrossprod(Xbar_re[,currentSample],Ubar_re[,currentSample]) + tcrossprod(Xbar_im[,currentSample],Ubar_im[,currentSample])) # E{xu^T}_re +~ (xbar_re %*% t(ubar_re) - xbar_im %*% t(-ubar_im)) ExuH_im <- (1-lambda) * ExuH_im + lambda * (tcrossprod(Xbar_im[,currentSample],Ubar_re[,currentSample]) - tcrossprod(Xbar_re[,currentSample],Ubar_im[,currentSample])) # E{xu^T}_im ~ (xbar_re %*% t(-ubar_im) + xbar_im %*% t(ubar_re)) } else if (complexModels) { ExuH_re <- (1-lambda) * ExuH_re + lambda * tcrossprod(Xbar_re[,currentSample],Ubar_re[,currentSample]) # E{xu^T}_re +~ xbar_re %*% t(ubar_re) ExuH_im <- (1-lambda) * ExuH_im - lambda * tcrossprod(Xbar_re[,currentSample],Ubar_im[,currentSample]) # E{xu^T}_im +~ xbar_re %*% t(-ubar_im) } else { ExuH_re <- (1-lambda) * ExuH_re + lambda * tcrossprod(Xbar_re[,currentSample],Ubar_re[,currentSample]) # E{xu^T}_re +~ xbar_re %*% t(ubar_re) } # E{xx^H}: if (complexSystem) { ExxH_re <- (1-lambda) * ExxH_re + lambda * (tcrossprod(Xbar_re[,currentSample]) + tcrossprod(Xbar_im[,currentSample])) # E{xx^T}_re +~ (xbar_re %*% t(xbar_re) - xbar_im %*% t(-xbar_im)) ExxH_im <- (1-lambda) * ExxH_im + lambda * (tcrossprod(Xbar_im[,currentSample],Xbar_re[,currentSample]) - tcrossprod(Xbar_re[,currentSample],Xbar_im[,currentSample])) # E{xx^T}_im +~ (xbar_re %*% t(-xbar_im) + xbar_im %*% t(xbar_re)) } else { ExxH_re <- (1-lambda) * ExxH_re + lambda * tcrossprod(Xbar_re[,currentSample]) # E{xx^T} +~ xbar_re %*% t(xbar_re) #invExxHold <- invExxH_re # experimenting iterative inversion by matrix inversion lemma #invExxH_re <- 1/(1-lambda) * ( invExxHold - (invExxHold %*% tcrossprod(Xbar_re[,currentSample]) %*% invExxHold) / c( t(Xbar_re[,currentSample]) %*% invExxHold %*% Xbar_re[,currentSample] + (1-lambda)/lambda) ) } } # Miscellaneous options: if (maskA) { # one way to use structural constraints for intelligent coupling: ExxH_re[row(ExxH_re)(imageH/2*imageW)) ] <- 0 #ExuH_re[ (row(ExuH_re)%%2 == 0) & (col(ExuH_re)<(imageH/2*imageW)) ] <- 0 if (intelligentCoupling) { # if intelligent coupling is assumed Q_re <- solve(ExxH_re) # set coupling matrix as inverse of activity covariance matrix } if (regulatedCoupling) { # if regulated coupling is assumed diag(Q_re) <- (1-lambda) * diag(Q_re) + lambda / pmax(diag(ExxH_re),0.01) # adapt each q_i towards 1/E{|x_i|^2} (or alternatively to 1/(3*E{|x_i|^2})), #diag(Q_im) <- (1-lambda) * diag(Q_im) + lambda / pmax(diag(ExxH_re),0.01) # adapt each q_i towards 1/E{|x_i|^2} (or alternatively to 1/(3*E{|x_i|^2})), #diag(R_re) <- (1-lambda) * diag(R_re) + lambda / pmax(diag(ExxH_re),0.01) # adapt each q_i towards 1/E{|x_i|^2} (or alternatively to 1/(3*E{|x_i|^2})), } # optimal coupling that maximizes E{|x_i|^2} (or E{|x_i||u_j|}) #Qp_re <- tcrossprod(sqrt(diag(Q_re))) #Q_re[row(Q_re) == col(Q_re)+1] <- 0.25*Qp_re[row(Q_re) == col(Q_re)+1] # add here possible interactions in q_ij #Q_re[row(Q_re) == col(Q_re)-1] <- 0.25*Qp_re[row(Q_re) == col(Q_re)-1] #Q_re[row(Q_re) == col(Q_re)+5] <- 0.25*Qp_re[row(Q_re) == col(Q_re)+5] # between neighbours on a 5x5 grid, for example #Q_re[row(Q_re) == col(Q_re)-5] <- 0.25*Qp_re[row(Q_re) == col(Q_re)-5] #Q_re <- tcrossprod(sqrt(diag(Q_re))) ##diag(R_re) <- 1/diag(Q_re) # response constants can be adapted too #diag(R_re) <- sqrt(diag(Q_re)) # response constants can be adapted too #varphi_re <- t(solve(solve(Q_re) + R_re %*% ExxH_re) %*% ExuH_re) # implicit adaptive filters as columns if (complexSystem) { phi_re <- t(Q_re %*% ExuH_re - Q_im %*% ExuH_im) phi_im <- -t(Q_re %*% ExuH_im + Q_im %*% ExuH_re) } else if (complexModels) { phi_re <- t(Q_re %*% ExuH_re) # implicit generative models as columns: compare to n eigen(Ubar%*%t(Ubar)/k)$vectors phi_im <- -t(Q_re %*% ExuH_im) # only in intelligent coupling case the first n models will be the same #varphi_im <- -t(solve(solve(Q_re) + R_re * ExxH_re) %*% ExuH_im) } else { phi_re <- t(Q_re %*% ExuH_re) # implicit generative models as columns: compare to n eigen(Ubar%*%t(Ubar)/k)$vectors phi_im <- 0 # or to n svd(Ubar,nu=n,nv=0)$u vectors } #U_im[1:n,] <- -Xbar_re # cross-correlations emerge if data is augmented with model activities # in addition, vectorizing successive samples together enables modeling of autocorrelations and power spectra if (showModels) { maxColor <- max(abs(range(phi_re,phi_im))) # calculate value range to center zero to white color #maxColor <- max(abs(range(phi_re*sqrt(varU),phi_im*sqrt(varU)))) # optionally reverse some of the preprocessing for (i in 1:n) { # draw all n models as images to observe emergence of structure values <- matrix(phi_re[,i],imageW,imageH) # models are in preprocessed data space #values <- matrix(basis %*% phi_re[,i]*sqrt(varU),imageW,imageH) # optionally reverse some of the preprocessing #values <- matrix(phi_re[,i]*sqrt(varU),imageW,imageH) # optionally reverse some of the preprocessing #values <- matrix(varphi_re[,i],imageW,imageH) # filters are in preprocessed data space image(values, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(bquote(paste(Re,group("{",phi[ .(i)],"}"))), side=1, col="black", cex=0.6) } if (complexModels) { values <- matrix(phi_im[,i],imageW,imageH) # models are in preprocessed data space #values <- matrix(phi_im[,i]*sqrt(varU),imageW,imageH) # optionally reverse some of the preprocessing #values <- matrix(varphi_im[,i],imageW,imageH) # filters are in preprocessed data space image(values, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(bquote(paste(Im,group("{",phi[ .(i)],"}"))), side=1, col="black", cex=0.6) } } } } if (showStructure) { if (complexModels) { phiHphi_re <- crossprod(phi_re) + crossprod(phi_im) # calculate inner products between models (t(phi_re) %*% phi_re - t(-phi_im) %*% phi_im) phiHphi_im <- crossprod(phi_re,phi_im) - crossprod(phi_im,phi_re) # calculate inner products between models (t(phi_re) %*% phi_im + t(-phi_im) %*% phi_re) } else { phiHphi_re <- crossprod(phi_re) # calculate inner products between models (t(phi_re)%*%phi_re) phiHphi_im <- 0 } maxColor <- max(abs(range(phiHphi_re,phiHphi_im))) # show model inner product matrix as an image image(phiHphi_re, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(expression(paste(Re,group("{",bold(phi)^H*bold(phi),"}"))), side=1, col="black", cex=0.6) } if (complexModels) { image(phiHphi_im, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(expression(paste(Im,group("{",bold(phi)^H*bold(phi),"}"))), side=1, col="black", cex=0.6) } } maxColor <- max(abs(range(ExxH_re,ExxH_im))) # show activity covariance matrix as an image image(t(ExxH_re), zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(expression(paste(Re,group("{",symbol(E):bold(bar(x))*bold(bar(x))^H,"}"))), side=1, col="black", cex=0.6) } if (complexSystem) { image(t(ExxH_im), zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(expression(paste(Im,group("{",symbol(E):bold(bar(x))*bold(bar(x))^H,"}"))), side=1, col="black", cex=0.6) } } } if (showCriterion) { p0 = 1 # study the minimization of cybernetic objective function if (p >= p0) { # but leave the first few values out to get a finer range # in case x_i is not cut to zero, these three objective functions are equivalent: # J[p-p0] <- 0.5 * sum( (Xbar_re %*% t(Xbar_re)) * (solve(Q_re) + R_re * ExxH_re) ) / k - sum( (Xbar_re %*% t(U)) * ExuH ) / k # J[p-p0] <- - 0.5 * sum( (Xbar_re %*% t(Xbar_re)) * (solve(Q_re) + R_re * ExxH_re) ) / k # J[p-p0] <- - 0.5 * sum( (Xbar_re %*% t(U)) * ExuH ) / k # also dampening of effective data could be studied # J[p-p0] <- sum( (Ubar %*% t(Ubar))^2 ) / k # but here we will just plot a simple Frobenius norm, not minimized, but indicative of adaptation if (!showStructure) { # calculate if not calculated already if (complexModels) { phiHphi_re <- crossprod(phi_re) + crossprod(phi_im) # calculate inner products between models (t(phi_re) %*% phi_re - t(-phi_im) %*% phi_im) phiHphi_im <- crossprod(phi_re,phi_im) - crossprod(phi_im,phi_re) # calculate inner products between models (t(phi_re) %*% phi_im + t(-phi_im) %*% phi_re) } else { phiHphi_re <- crossprod(phi_re) # calculate inner products between models (t(phi_re)%*%phi_re) phiHphi_im <- 0 } } #J[p-p0] <- sum(phiHphi_re^2+phiHphi_im^2) # square of Frobenius norm, or plot only the diagonal similar to sum(|ExuH|^2) or some custom functional #J[p-p0] <- sum(phiHphi_re^2+phiHphi_im^2) - 2*sum(diag(phiHphi_re)^2) # modified to reverse the cost of diagonal Jcur <- sum(phiHphi_re) - 2*sum(diag(phiHphi_re)) # sum of entries, modified to reverse the cost of diagonal J[p-p0] <- Jcur #(1-abs((Jcur-Jprev)/Jprev))^2 #log(1/abs(Jcur - Jprev))/40 #Jprev <- Jcur plot(J, type="l", bty="l", yaxt="n") # or plot only the latest values or some logarithmic transform to see more nuances } } if (showEstimate) { someSample = (p-1)%%k + 1 # index of a sample, here iteration number p limited to number of samples k maxColor <- maxU * 1.5 # set value range large enough for data, center zero to white color values <- matrix(U_re[,someSample]*sqrt(varU)+meanU_re,imageW,imageH) # draw one sample as image, reverse preprocessing #values <- matrix(U_re[,someSample],imageW,imageH) # alternatively draw one sample as image as is, update maxColor image(values, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(bquote(paste(Re,group("{",bold(u)*(.(someSample)),"}"))), side=1, col="black", cex=0.6) } if (complexModels) { values <- matrix(U_im[,someSample]*sqrt(varU)+meanU_im,imageW,imageH) # draw one sample as image, reverse preprocessing #values <- matrix(U_im[,someSample],imageW,imageH) # alternatively draw one sample as image as is image(values, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(bquote(paste(Im,group("{",bold(u)*(.(someSample)),"}"))), side=1, col="black", cex=0.6) } } values <- matrix((U_re[,someSample]-Ubar_re[,someSample])*sqrt(varU)+meanU_re,imageW,imageH) # draw estimate as image #values <- matrix(phi_re %*% R_re %*% Xbar_re[,someSample]*sqrt(varU)+meanU_re,imageW,imageH) # draw estimate as image #values <- matrix(phi_re %*% R_re %*% Xbar_re[,someSample],imageW,imageH) # draw implicit estimate as image #values <- (values > 0.5) # optionally threshold estimates back to binary data image(values, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(bquote(paste(Re,group("{",bold(hat(u))*(.(someSample)),"}"))), side=1, col="black", cex=0.6) } if (complexModels) { values <- matrix((U_im[,someSample]-Ubar_im[,someSample])*sqrt(varU)+meanU_im,imageW,imageH) # draw estimate as image #values <- matrix(phi_im %*% R_re %*% Xbar_re[,someSample]*sqrt(varU)+meanU_im,imageW,imageH) # draw estimate as image #values <- matrix(phi_im %*% R_re %*% Xbar_re[,someSample],imageW,imageH) # draw implicit estimate as image #values <- (values > 0.5) # optionally threshold estimates back to binary data image(values, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(bquote(paste(Im,group("{",bold(hat(u))*(.(someSample)),"}"))), side=1, col="black", cex=0.6) } } } if (showBalance) { maxColor <- maxUbar * 1.5 # set value range large enough for data, center zero to white color #values <- matrix(Ubar_re[,someSample]*sqrt(varU)+meanU_re,imageW,imageH) # draw balanced data as image, reverse preprocessing values <- matrix(Ubar_re[,someSample],imageW,imageH) # alternatively draw one balance as image as is, update maxColor image(values, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(bquote(paste(Re,group("{",bold(bar(u))*(.(someSample)),"}"))), side=1, col="black", cex=0.6) } if (complexModels) { #values <- matrix(Ubar_im[,someSample]*sqrt(varU)+meanU_im,imageW,imageH) # draw balanced data as image, reverse preprocessing values <- matrix(Ubar_im[,someSample],imageW,imageH) # alternatively draw one balance as image as is image(values, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(bquote(paste(Im,group("{",bold(bar(u))*(.(someSample)),"}"))), side=1, col="black", cex=0.6) } } values <- matrix(c(Xbar_re[,someSample],rep(0,gridSize^2-n)),gridSize,gridSize) # draw one activity as image, update maxColor maxColor <- max(abs(range(values))) image(values, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(bquote(paste(Re,group("{",bold(bar(x))*(.(someSample)),"}"))), side=1, col="black", cex=0.6) } if (complexSystem) { values <- matrix(c(Xbar_im[,someSample],rep(0,gridSize^2-n)),gridSize,gridSize) # draw one activity as image, update maxColor maxColor <- max(abs(range(values))) image(values, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) if (showLabels) { mtext(bquote(paste(Im,group("{",bold(bar(x))*(.(someSample)),"}"))), side=1, col="black", cex=0.6) } } } # optionally save results as numbered images to current working directory #savePlot(filename = paste("image",formatC(p,width=5,format="d",flag="0"), sep=""),type="png") } # after convergence, study image(Xbar_re), image(U), image(Ubar), # image(phi), image(phiHphi_re) vs image(Q_re%*%ExxH_re), image(Uest) as R_re*phi%*%Xbar_re # (or (R_re+1)*phi%*%Xbar_re to double the magnitude and allow R_re=0), # their ranges, diagonals, covariances and eigenvalues, etc. # remember to cut negative and imaginary Xbar to zero and # to reverse preprocessing to use estimates in the original data space