31:T67c,md` ## Goal Walk through each of the steps needed to perform t-SNE. We will execute each step explicitely in plain javascript, no libraries allowed. At the end there will be a bunch of sliders that allow you to play with the algorithm and watch it attempt to visualize some simulated high-dimensional data. --- ## Calculate pairwise distances. The first step of t-SNE is to calculate the pairwise distances between all of the points in your data. Luckily we only need to do this a single time. This is because t-sne only cares about these distances and they don't change ever change. This only-once property is great in that after you have done the distance calculations, the actual data doesn't matter. So really the runtime of the algorithm is bounded not by the dimensionality of your input data but by how many samples you have. The nieve way of doing this pairwise calculation would be a simple outer loop scanning over each observation in turn and then within that outer loop an inner loop that calculates the distance between it and every other point. This works but is wasteful. Let's think about these distances. The distance between point \`i\` and point \`j\` is the same as the distance between \`j\` and point \`i\` so really we should only calculate a single one, in addition we don't really need to check between \`i\` and \`i\` as we know it's zero. To see how much time this can save us we can use the formula for calculating all _permutations_(order matters) of two numbers and all _combinations_(order doesn't matter) of two numbers. With ${n} total datapoints or observations, we can calculate the total number of permutations: `32:T892,function beta_binary_search(p_row, i, pw_dists, entropy_tol = 1e-3, max_attempts = 40){ let beta_min = -Infinity, // beta upper bound beta_max = Infinity, // beta lower bound beta = 1, // initial beta value searching_for_beta = true, attempts = 0; // number of iterations through while looking while(searching_for_beta){ // binary search to find precision beta so our output entropy is what we want. // we need to zero out our sum of probabilities and entropy at each iteration let sum_p = 0; let entropy = 0; // initial entropy for this loop attempts++ // increment out attempts counter // first we calculate the vector of probabilties for pi_{all other points} for(let j=0; j 1e-7){ entropy -= p_row[j]*log(p_row[j]) } } // check where our entropy lies in context to the target if(entropy > target_entropy){ // distribution is too even. beta_min = beta; // move our min bound up // set new beta based on our bounds. beta = beta_max === Infinity ? beta*2 : (beta+beta_max)/2; } else { // distribution is too clumped. beta_max = beta; // move our max bound down // set new beta based on our bounds. beta = beta_min === -Infinity ? beta/2 : (beta+beta_min)/2; } const too_many_attempts = attempts > max_attempts; const close_to_target = Math.abs(entropy - target_entropy) < entropy_tol; searching_for_beta = !(too_many_attempts || close_to_target) } return p_row }33:T79c,Y_costs = { let iters = 0 const costs = []; // initial random y positions along with a vector to store the last y gradient step and the gains. const Y_new = initialize_points({n: num_points, dim: map_dims}), Y_step = initialize_points({n:num_points, dim: map_dims, fill:0}), Y_gains = initialize_points({n: num_points, dim: map_dims, fill: 1}); let average_y = make_array(map_dims); // calculate the cost and gradient for current y locations while(iters < num_itterations){ iters++ average_y = average_y.map(() => 0) // zero out y average for new iteration let {cost, grad} = calc_cost_grad(Y_new, iters), momentum = momentums[iters > momentum_switch ? 1: 0], grad_id, step_id, gain_id, gain_id_new, step_id_new; costs.push(cost); for(let i=0; i { const low = -cluster_spread/2; const high = cluster_spread/2; const point_gen = d3.randomUniform(low, high); return () => [...(new Array(ndim))].map(point_gen); } const gen_centroids = (n, ndim) => { const center_maker = make_centroid(ndim); return [...(new Array(n))].map(center_maker) } const centroids = gen_centroids(num_clusters, num_dimensions) const sample_from_centroid = (nsamples, centroid, centroid_id) => { const rnorm = d3.randomNormal(0, random_noise_sd); const get_sample = () => ({group:centroid_id, location : centroid.map(d => d + rnorm())}); return [...(new Array(nsamples))].map(get_sample); } function gen_data_from_centroids(centroids, nsamples){ return centroids.reduce( (data, d, i) => [...data, ...sample_from_centroid(nsamples, d, i)], []) .map((d,i) => (Object.assign({id: i},d,{}))) // add in node id in addition to the centroid one. } return gen_data_from_centroids(centroids, num_samples) }35:[] 2f:{"title":"t-SNE in plain javascript","cells":[{"id":0,"value":"md`# t-SNE in plain javascript`","pinned":false,"hidden":false,"mode":"ojs"},{"id":54,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":835,"value":"viewof n = slider2('Set number of observations in your dataset', 'n', 0, 100, 50)","pinned":false,"hidden":false,"mode":"ojs"},{"id":785,"value":"tex`\n\\begin{aligned}\nn^2 = ${n}^2 = \\boxed{${squared(n)} \\text{ pairs}}\n\\end{aligned}`\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":789,"value":"md`\nNow if we're efficient about it and only look at the combinations of the data we can cut this down. The combination equation is a bit more complicated...\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":798,"value":"tex`\n\\begin{aligned}\n \\frac{n!}{2!(n - 2)!} =& \\frac{${n}!}{2\\cdot${n - 2}!} \\\\\n& \\frac{ ... ${n - 3} \\cdot ${n - 2} \\cdot ${n - 1} \\cdot ${n}}{2\\cdot(...${n - 3} \\cdot ${n - 2})} \\\\\n& \\frac{${n - 1} \\cdot ${n}}{2} \\\\\n& \\boxed{${n*(n - 1)/2}\\text{ pairs}}\n\\end{aligned}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":797,"value":"md`So we can see that by taking advantage of this combination vs permutation distinction we are saving ourselves __${squared(n) - (n*(n - 1)/2)}__ total distance calculations. Kinda nice, especially when the data gets larger. \n\nNow that we've gotten some superfluous math out of the way let's actually implement this in code. A couple things to note:\n\n1. The function \\`calc_distance\\` simply takes two arrays of some dimension and calculates their euclidean distance(sometimes called the _L2_ norm) from eachother. The actual code to implement it is at the bottom of this notebook if you're curious.\n2. We are filling an array that has every simple permutation of the data in it, but if you note the second for loop's bounds we are only calculating the combinations. We simple fill in two spots in our distances array for each distance. This is because we later need the array to be the size of all the distances due to the quirk of non-symmetry in the conditional probability step. (We'll get to it, don't worry.)\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":62,"value":"pairwise_distances = {\n let x_i, x_j, distance;\n const pw_dists = new Array(num_points*num_points);\n for(let i = 0; i < num_points; i++){\n \tx_i = data[i].location;\n for(let j = i+1; j < num_points; j++){\n x_j = data[j].location;\n distance = calc_distance(x_i, x_j);\n pw_dists[i*num_points+j] = {source: i, target: j, distance};\n pw_dists[j*num_points+i] = {source: j, target: i, distance};\n }\n }\n return pw_dists\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":59,"value":"md`Now we need to calculate the probabilities associated with the high dimensional spaces`","pinned":false,"hidden":false,"mode":"ojs"},{"id":114,"value":"target_entropy = log(perplexity) // entropy level we're shooting for given perplexity","pinned":true,"hidden":false,"mode":"ojs"},{"id":237,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":273,"value":"md`Now we calculate the non-symmetric probabilities. Aka probability of \\`i\\` given \\`j\\`.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":108,"value":"P_i_given_j = {\n const P_i_given_j = make_array(num_points*num_points); \n const p_row = make_array(num_points);\n \n for(let i = 0; i < num_points; i++){\n // run binary search to find the best beta. \n let new_p_row = beta_binary_search(p_row, i, pairwise_distances)\n for(let j=0;j squared(Y_curr[i][dim] - Y_curr[j][dim]));\n \n q_ij = 1/(1 + dim_sum); // set the q value for this combination\n Q[i*num_points+j] = q_ij; // assign it to the big Q array\n Q[j*num_points+i] = q_ij;\n sum_q += 2*q_ij; // accumulate for later normalization\n }\n }\n // normalize Q values\n return Q.map(q => max(q/sum_q, 1e-99));\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":355,"value":"function calc_cost_grad(Y, iters){\n let cost_acc = 0,\n pre_mult, // this is poorly named;\n grad_sum,\n early_exag = iters < early_exag_itts? early_exag_amount: 1;\n\n const Q = find_Q(Y),\n \t\tgrad_acc = make_array(num_points);\n \n for(let i=0; i0); // initialize array of zeros for gradient accumulation\n for(let j=0; j grad + (pre_mult * (Y[i][dim] - Y[j][dim])))\n }\n \n grad_acc[i] = grad_sum;\n }\n return {cost: cost_acc, grad: grad_acc}\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":631,"value":"calc_cost_grad(initialize_points({n: num_points, dim: map_dims}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":410,"value":"md`Now that we have the cost and gradient at this step we can recalculate our y mapping by moving in the direction of the gradient until we lower our cost as much as possible.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":412,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":731,"value":"Y = Y_costs.Y","pinned":false,"hidden":false,"mode":"ojs"},{"id":733,"value":"costs = Y_costs.costs","pinned":false,"hidden":false,"mode":"ojs"},{"id":581,"value":"{\n const svg = d3.select(DOM.svg(width, height));\n const dim1_domain = d3.extent(Y, (y) => y[0]);\n const dim2_domain = d3.extent(Y, (y) => y[1]);\n \n const dim1 = d3.scaleLinear().domain(dim1_domain).range([margin.left, width - margin.right]);\n const dim2 = d3.scaleLinear().domain(dim2_domain).range([margin.top, height - margin.bottom]);\n const group_color = d3.scaleOrdinal()\n \t.domain(make_array(num_clusters).map((d,i)=>i))\n .range(d3_color.schemeCategory10);\n \n const vis_data = data.map((d,i) => ({...d, map: Y[i]}));\n \n svg.append(\"g\")\n .attr(\"stroke\", \"#000\")\n .attr(\"stroke-opacity\", 0.2)\n .selectAll(\"circle\")\n .data(vis_data)\n .enter().append(\"circle\")\n .attr(\"cx\", d => dim1(d.map[0]))\n .attr(\"cy\", d => dim2(d.map[1]))\n .attr(\"fill\", d => group_color(d.group))\n .attr(\"r\", 5);\n \n return svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":729,"value":"{\n const svg = d3.select(DOM.svg(width, height));\n \n const x = d3.scaleLinear().domain([0, costs.length]).range([margin.left, width - margin.right]);\n const y = d3.scaleLinear().domain(d3.extent(costs)).range([height - margin.bottom, margin.top]);\n\n svg.append(\"g\")\n .attr(\"stroke\", \"#000\")\n .attr(\"stroke-opacity\", 0.2)\n .selectAll(\"circle\")\n .data(costs)\n .enter().append(\"circle\")\n .attr(\"cx\", (d,i) => x(i))\n .attr(\"cy\", y)\n .attr(\"fill\", 'orangered')\n .attr(\"r\", 2);\n \n return svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":725,"value":"viewof num_itterations = slider({min: 1, max: 10000, value: 1000})","pinned":false,"hidden":false,"mode":"ojs"},{"id":780,"value":"viewof perplexity = slider({min: 1, max: 70, value: 30})","pinned":false,"hidden":false,"mode":"ojs"},{"id":770,"value":"viewof learning_rate = slider({min: 0, max: 20, value: 10})","pinned":false,"hidden":false,"mode":"ojs"},{"id":574,"value":"viewof early_exag_amount = slider({min: 1, max: 10, value: 4})","pinned":false,"hidden":false,"mode":"ojs"},{"id":567,"value":"viewof early_exag_itts = slider({min: 0, max: num_itterations, value: 250})","pinned":false,"hidden":false,"mode":"ojs"},{"id":492,"value":"viewof gain_scales_same = slider({min: 0, max: 1, value: 0.8})","pinned":false,"hidden":false,"mode":"ojs"},{"id":451,"value":"viewof gain_scales_diff = slider({min: 0, max: 1, value: 0.2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":458,"value":"viewof gain_clamp = slider({min: 0, max: 0.1, value: 0.01})","pinned":false,"hidden":false,"mode":"ojs"},{"id":502,"value":"viewof momentum_switch = slider({min: 0, max: 1000, value: 250})","pinned":false,"hidden":false,"mode":"ojs"},{"id":466,"value":"momentums = [0.5, 0.8]","pinned":false,"hidden":false,"mode":"ojs"},{"id":483,"value":"map_dims = 2","pinned":false,"hidden":false,"mode":"ojs"},{"id":428,"value":"md`Logic for generating data is here, It's not interesting so I haven't pinned it.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":707,"value":"num_dimensions = 10","pinned":false,"hidden":false,"mode":"ojs"},{"id":711,"value":"num_clusters = 4","pinned":false,"hidden":false,"mode":"ojs"},{"id":714,"value":"num_samples = 15","pinned":false,"hidden":false,"mode":"ojs"},{"id":864,"value":"num_points = data.length","pinned":false,"hidden":false,"mode":"ojs"},{"id":716,"value":"cluster_spread = 15","pinned":false,"hidden":false,"mode":"ojs"},{"id":742,"value":"viewof random_noise_sd = html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":718,"value":"random_noise_sd","pinned":false,"hidden":false,"mode":"ojs"},{"id":11,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":538,"value":"function calc_distance(x1, x2){\n // calculates the euclidean distance between two vectors. Aka the L2 Norm\n const sum_of_square_diffs = x1.map(\n (x1_i, i) => Math.pow(x1_i - x2[i], 2) // get squared difference each dimension\n ) .reduce((sum, d) => sum + d, 0); // sum the squared differences to return a single number\n \n return Math.sqrt(sum_of_square_diffs)\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":140,"value":"exp = Math.exp","pinned":false,"hidden":false,"mode":"ojs"},{"id":341,"value":"squared = (val) => Math.pow(val, 2)","pinned":false,"hidden":false,"mode":"ojs"},{"id":352,"value":"max = Math.max","pinned":false,"hidden":false,"mode":"ojs"},{"id":147,"value":"log = Math.log","pinned":false,"hidden":false,"mode":"ojs"},{"id":439,"value":"sign = (x) => x > 0 ? 1: x < 0 ? -1: 0","pinned":false,"hidden":false,"mode":"ojs"},{"id":551,"value":"function initialize_points({n, dim = 2, fill}){\n const fill_function = typeof fill == 'undefined' ?\n d3.randomUniform(-1,1) :\n \t\t() => fill;\n const rand_point = () => make_array(dim).map(fill_function);\n return make_array(n).map(rand_point)\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":327,"value":"make_array = (n) => [...(new Array(n))]","pinned":false,"hidden":false,"mode":"ojs"},{"id":703,"value":"margin = ({top: 20, right: 30, bottom: 30, left: 40})","pinned":false,"hidden":false,"mode":"ojs"},{"id":583,"value":"height = 300","pinned":false,"hidden":false,"mode":"ojs"},{"id":721,"value":"d3_color = require('d3-scale-chromatic')","pinned":false,"hidden":false,"mode":"ojs"},{"id":849,"value":"function slider2(msg, var_name, min, max, start){\n const div = html`\n\t${msg}\n \n ${var_name} = \n `;\n const range = div.querySelector(\"[type=range]\");\n const number = div.querySelector(\"[type=number]\");\n div.value = range.value = number.value = start;\n range.addEventListener(\"input\", () => number.value = div.value = range.valueAsNumber);\n number.addEventListener(\"input\", () => range.value = div.value = number.valueAsNumber);\n return div;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":751,"value":"import {slider} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":47,"value":"d3 = require('d3')","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}