# Raw procedure for learning emergent models by computing cybernetic equilibria. # Run in R statistical computing environment with 'source("sca.r")' in proper directory. # With intelligent coupling Q = E{xx^T}^-1 between systems x_i, and with zero response/consuming constants R, # the models phi_i will 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 coupling q_i = constant # by consuming the data for activity, R > 0, which separates the models. In addition, constraint x_i > 0 # elicits sparse components. With careful conjugations, also complex values can be implemented. # The models summarize steady states of rapid gradient descents of cybernetic objective function # J(x,u) = 1/2 x^T A x - x^T B u, i.e. dx/dt = - G A x + G B u, i.e. xbar = A^-1 B u = phi^T ubar, # where A = (Q^-1 + R E{xbarxbar^T}) and B = E{xbarubar^T}, # 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 <- 15 # number of models/features/systems/attractors x_i lambda <- 0.1 # near zero: inverse of adaptation time scale in sample time; reduce if problems 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 centerDimensions <- TRUE # center data (recommended for numerical stability) normalizeDimensions <- FALSE # normalize data variances to unity (recommended for unknown data) initModels <- TRUE # initialize model matrices useRays <- TRUE # enhance sparsity by cutting negative inner products to zero maskA <- FALSE # elicit actual principal components by structural constraints q_i <- 4 # coupling constants' initial values intelligentCoupling = FALSE # Q = E{xx^T}^-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 r_i <- 1 # response constants between 0 (decoupled) and 1 (consuming data for activity) showModels = TRUE # show models as images showStructure = TRUE # show higher-order structure as images, matrices phi^T phi and E{xx^T} for example showCriterion = TRUE # plot some objective function of iteration, showing convergence for example showEstimate = TRUE # show some data estimates by current models if (loadData) { m <- 256 # number of data dimensions u_j k <- 500 # number of data samples u U <- matrix(scan("digits500.txt"), m, k) # load data samples to columns of matrix U # try library 'matlab' or 'R.matlab' for .mat data files side <- 16 # image width and height for this special data if (centerSamples) { # optionally center sample vectors sampleMeans <- colMeans(U) # calculate sample means U <- t(t(U)-sampleMeans) # center samples by substracting the means } if (normalizeSamples) { # optionally normalize sample vectors to unit length sampleSqMeans <- colMeans(U^2) # calculate means of sample squares (std.dev^2, RMS^2) sampleSqMeans[sampleSqMeans==0] <- 1 # prevent division by zero U <- t(t(U)/sqrt(sampleSqMeans)) # normalize sample vectors to unit length } maxU <- max(abs(range(U))) # store maximum value for drawing if (centerDimensions) { # for numerical stability it is better to center the data meanU <- rowMeans(U) # calculate data dimension averages to set affine subspaces U <- U - meanU # center the data by substracting the mean (store for reversal) } else { meanU <- 0 } if (normalizeDimensions) { # if nothing more specific is known about the data distribution varU <- rowMeans(U^2) # calculate data dimension variances #varU <- mean(varU) # optionally treat several dimensions as similar varU[varU==0] <- 1 # prevent division by zero U <- U/sqrt(varU) # normalize data by dividing by standard deviations } else { varU <- 1 } } # for natural positive data consider instead logarithms and # geometric means of data, for example dividing deviations # from means by those means. # for time series study vectorizing successive samples together if(initModels) { ExuT <- 0.1*matrix(runif(n*m)-0.5,n,m) # initialize matrix E{xu^T} with small random numbers (1/sqrt(m), for example) ExxT <- 0.1*matrix(runif(n*n),n,n) # initialize matrix E{xx^T} with small positive random numbers (or just a diagonal) Q <- q_i*diag(n) # initialize matrix Q with equal coupling constants in diagonal if (intelligentCoupling) { # if intelligent coupling is assumed Q <- solve(ExxT) # initialize coupling matrix as inverse of activity covariance matrix } phi <- t(Q %*% ExuT) # initialize models as columns of matrix phi #phiTphi <- crossprod(phi) # calculate inner products between models (t(phi)%*%phi) R <- r_i*diag(n) # initialize matrix R with equal response constants in diagonal } # 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 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 for (p in 1:1000) { # iterate for some time or when Esc pressed par(mfrow=c(5,5), mar=c(1,1,1,1)) # create/clear a graphics window with 5x5=25 tiled image blocks Xbar <- solve(solve(Q) + R %*% ExxT) %*% ExuT %*% U # calculate steady state for each u ##Xbar <- solve(diag(n) + R %*% phiTphi) %*% t(phi) %*% U # equivalent alternative formulation if (useRays) { # one way to elicit sparse solutions: Xbar[Xbar<0] <- 0 # cut negative inner products to zero # proper nonlinearity should be applied recursively with active models ##Xbar <- t(phi) %*% U - R %*% phiTphi %*% Xbar # but not possible in a matrix form ##Xbar[Xbar<0] <- 0 # one could study applying these approximations a few times } Ubar <- U - phi %*% R %*% Xbar # calculate the effective data by substracting the effect of # steady states: the estimate that is implicit in consuming responses # update system/model matrices ExuT <- (1-lambda) * ExuT + lambda * tcrossprod(Xbar,Ubar)/k # E{xu^T} ~ Xbar %*% t(Ubar) / k ExxT <- (1-lambda) * ExxT + lambda * tcrossprod(Xbar)/k # E{xx^T} ~ Xbar %*% t(Xbar) / k if (maskA) { # one way to use structural constraints for intelligent coupling: ExxT[row(ExxT)= p0) { # but leave the first few values out to get a finer range for y-axis # in case x_i is not cut to zero, these three objective functions are equivalent: # J[p-p0] <- 0.5 * sum( (Xbar %*% t(Xbar)) * (solve(Q) + R %*% ExxT) ) / k - sum( (Xbar %*% t(U)) * ExuT ) / k # J[p-p0] <- - 0.5 * sum( (Xbar %*% t(Xbar)) * (solve(Q) + R %*% ExxT) ) / k # J[p-p0] <- - 0.5 * sum( (Xbar %*% t(U)) * ExuT ) / 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) { phiTphi <- crossprod(phi) } # calculate if not calculated already J[p-p0] <- sum(phiTphi^2) # or plot only the diagonal similar to sum(ExuT^2), or some custom functional plot(J, type="l", bty="n", yaxt="n") # or plot only 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[,someSample]*sqrt(varU)+meanU,side,side) # draw one sample as image, reverse preprocessing image(values, zlim=c(-maxColor,maxColor), col=bwr, axes=FALSE, ylim=c(1.05,-0.05)) values <- matrix(phi %*% (R+diag(n)) %*% Xbar[,someSample]*sqrt(varU)+meanU,side,side) # draw 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)) } } # after convergence, study Xbar, U, phi, Uest as phi%*%R%*%Xbar # (or phi%*%(I+R)%*%Xbar to double the magnitude and allow R=0), # their covariances and eigenvalues, etc. # remember to cut negative Xbar to zero and to reverse preprocessing # to use estimates in the original data space