Multidimensional Scaling (MDS) is a linear embedding method used when we only know the pairwise distances between data points. For linear systems, MDS works well with as little as 10 points and the system is 2 dimensional. For this program, you choose the window size or embedding for the system and MDS identifies what the true embedding should be. The following program uses Logistic map data when *r = 3.2* but you can change *r* to see how well MDS works as the map bifurcates and eventually becomes chaotic and nonlinear at 4. MDS breaks down at this point and the embedding defaults to just 1 over the embedding you choose.

# Requires Numpy and Scipy libraries from numpy import * import scipy.linalg import random r=3.2 # parameter on Logistic map n = 10 # the number of timeseries points to use d = 40 # your embedding choice # Create a time series def timeseries(number,throw, d): xList = [] dataX0 = [] t = 0 while t < d: xList.append(.1) t += 1 while t <= number+throw+d: x = r * xList[t-1] * (1 - xList[t-1]) # Logistic Map xList.append(x) if t > throw: y = xList[t-1] dataX0.append(x) t = t + 1 return dataX0 # Construct an n x n centering matrix # The form is P = I - (1/n) U where U is a matrix of all ones def centering_matrix(n): P = eye(n) - 1/float(n) * ones((n,n)) return P def create_data(number, d): dataX0 = timeseries(number, 500, d) # Create a sliding window of size d data = [] i = d while i < len(dataX0): data.append([]) j = i - d while j < i: data[len(data)-1].append(dataX0[j]) j = j + 1 i = i + 1 return data def main(): print "n =",n, ": number of timeseries points used" print "d =",d, ": embedding" P = centering_matrix(n) # Create the data with the number of dimensions data = create_data(n, d) X = pairwise_distances(data) A = -1/2.0 * P * X * P # Calculate the eigenvalues/vectors [vals, vectors] = scipy.linalg.eig(A) # Sort the values vals = sort(vals) embedding = 0 for x in vals: if x > 10**(-10): embedding += 1 print "mds embedding =",embedding # Compute the pairwise distance between vector x and y def metric(x1, y1): d = 2 summ = [] i = 0 while i < len(x1): # in this case use euclidean distance summ.append((x1[i] - y1[i])**d) i = i + 1 return sum(summ) ** (1 / float(d)) # Return a matrix of pairwise distances def pairwise_distances(data): distances = [] i = 0 while i < len(data): distances.append([]) x1 = data[i] j = 0 while j < len(data): y1 = data[j] distances[i].append(metric(x1, y1)**2) j = j + 1 i = i + 1 return distances main()