2f:T110b,md` # Goal This notebook attempts to be a guide that walks you through every step of the t-SNE algorithm so that when you use it you understand what it is doing entirely. It will introduce each concept/formula associated with t-SNE and then implement it in javascript. Everything is totally editable so feel free to tinker with the functions and see if you break or improve something! Each step will be written in plain plain javascript, no libraries allowed (although I do make visualizations using D3). 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. Not a javascript person? I have attempted to make the code as universal as possible (e.g. favoring for loops over some javascript specific vectorized functions) so the general ideas will hopefully transfer! _A large source of help in coding this came from the excellent [tSNEjs](https://github.com/karpathy/tsnejs) library made by [Andrej Karpathy.](https://twitter.com/karpathy)_ *Once you're done with this I highly recomend exploring some of the pitfalls with the excellent interactive article [How to Use t_SNE Effectively.](https://distill.pub/2016/misread-tsne/)* ## What t-SNE is t-SNE stands for "[T-Distributed Stochastic Neighbor Embedding](https://lvdmaaten.github.io/tsne/)" and is a method for taking data of a high enough dimension that it can't be visualized in its raw form and maps it to a lower dimension (usually 2 or 3 dimensions). This has wide ranging applications wherever discovering structure in data is valuable. For example: [looking at the expressed proteins in a given cell population to attempt to find different cell types](https://tinyurl.com/ybgs9fv9). ## What t-SNE _isn't_ It is not a dimensionality reduction technique like [Principle Components Analysis](https://en.wikipedia.org/wiki/Principal_component_analysis). This means you can't train t-SNE on data and then feed new data into it to get predictions about where they are in the mapping. One major implication of this is it can't be used to say, map data to a lower dimension and then use a clustering algorithm to classify. There are however some attempts at taking the benefits of t-SNE and also having it be a dimensionality reduction technique such as [parametric t-SNE](https://lvdmaaten.github.io/publications/papers/AISTATS_2009.pdf). --- # The t-SNE steps Roughly we can think about t-SNE in 4 steps: 1. Calculate pairwise distances between all points in their high dimension. 2. Convert those pairwise distances to pairwise-probabilities. 3. Map each point to a lower dimensional space 4. Iterate this mapping. until pairwise probabilities in the lower-dimensional mapping are close to the pairwise probabilities in the high dimension. Okay, let's look at these steps in detail. ## Calculate pairwise distances. The first step of t-SNE is to calculate the pairwise distances between all of the points in your data in their original high-dimensional form. One nice think about the fact that t-SNE only cares about pairwise distances is that we only need to calculate these once and then never touch them again. So really the runtime of the algorithm is bounded not by the dimensionality of your input data but by how many samples you have. After calculating pairwise distances t-SNE would run as fast on ${tex`64 \times 64`} dimensional MNIST images as it would on our simulated ${num_dimensions} dimensional data. The naieve 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 ${tex`X_i`}$ and point ${tex`X_j`}$ is the same as the distance between ${tex`X_j`} and point ${tex`X_i`} so really we should only calculate a single one, in addition we don't really need to check between ${tex`X_i`} and ${tex`X_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: `30:T484,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. Now that we've gotten some superfluous math out of the way, let's actually implement this in code. A couple things to note: 1. We have a few helper functions that will help clean up our code a bit. These are \`calc_distance\` which calculates the euclidean distance (sometimes called the _L2_ norm) between two arrays, \`mat_int\` which helps us translate from a 2-d matrix position (row i, column j) to a 1-d javascript array, and lastly \`squared\` which does what it sounds like, squares its supplied value. 2. 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.) `31:T5cf,md` If you expand the returned array from above you will see we have a series of objects containing \`source\`: the index in our original data corresponding to the first point, \`target\`: the index for our second point, and \`distance\`: the euclidean distance between those two points. ## Pairwise probability calculations Now that we have the distances we can get to the fun stuff. We need to calculate the conditional probabilities between points based on their distance from eachother in the original high dimensional space. ### Going from distance to probability We can think about this as taking a given point \`i\` and then calculating the distance to every other point from it. We then place a normal or gaussian distribution on top of the point (or simply zero distance) and then use that to calculate the 'probability' of every other point based on their distance from the original. This is a bit to wrap your head around so let's look at it in the case of just two points. In the following plot we are showing the calculation of the conditional probability of point i given a gaussian centered around point j. Play around with the position of point i and the standard deviation of the distribution to see the relationship between euclidean distances between points and their conditional probabilies.`32:T663,{ const pdf = linspace(-5, 5, 350) .map(x => ({x, px:stats.dnorm(x, 0, orig_variance)})), x = d3.scaleLinear() .domain(d3.extent(pdf, d => d.x)) .range([margin.left, width - margin.right]), y = d3.scaleLinear() .domain(d3.extent(pdf, d => d.px)) .range([height - margin.bottom, margin.top]), area = d3.area() .x(d => x(d.x)) .y1(y(0)) .y0(d => y(d.px)), p_ij = stats.dnorm(new_point, 0, orig_variance); const svg = d3.select(DOM.svg(width, height)); svg.append("path") .datum(pdf) .attr("stroke", "darkgrey") .attr("stroke-width", 1.5) .attr("stroke-linejoin", "round") .attr("stroke-linecap", "round") .attr("fill", "steelblue") .attr("fill-opacity", 0.2) .attr("d", area); svg.append('circle') .attr('cx', x(0)) .attr('cy', y(0)) .attr('r', 5) .attr('fill', 'steelblue'); svg.append('circle') .attr('cx', x(new_point)) .attr('cy', y(0)) .attr('r', 5) .attr('fill', 'orangered'); svg.append("line") .attr("fill", "none") .attr("stroke", "orangered") .attr("stroke-width", 1.5) .attr("stroke-linejoin", "round") .attr("stroke-linecap", "round") .attr("x1", x(new_point)) .attr("x2", x(new_point)) .attr("y1", y(p_ij)) .attr("y2", y(0)) svg.append("text") .html( `P(i|j) = ${p_ij.toPrecision(3)}` ) .attr('font-size', 20) .attr("x", x(new_point) + 3) .attr("y", y(p_ij) + 10) return svg.node(); }33:T5ea,md` ## Perplexity You may notice that in the above example we can tune the variance, or width, of the distribution that we center over a point to calculate it's respective conditional probabilities to other points. The perplexity hyper-parameter of t-SNE helps control this. ### Why is it perplexity and not just variance? It's tempting to say that in the final algorithm we just set the same variance for every point's assigned distribution. Unfortunately that doesn't work very well. The reason is we may have areas in the space that are very dense with points and others that are very sparse. Since we care about structure at multiple scales treating data in both of these regions the same will not give us good results. The way that t-SNE deals with this is by scanning over possible standard deviations for every point's distribution trying to fit to a parameter called 'perplexity'. ### What is it? Perplexity a measure related to the [entropy](https://en.wikipedia.org/wiki/Entropy) (or dispersed-ness) of the system of points. A valuable way of thinking about perplexity is to think of it as the effective number of neighbors for each point. Aka how many points that the distribution over point \`i\` capture with a non-trivial probability. The equation is relatively simple. (_The calculations for ${tex`p_{i|j}`} are coming up next, so don't worry about them yet_)... ${tex.block`\log(\text{perplexity for point } i) = (\text{entropy for point } i) = \sum_{j \ne i} -p_{j|i}\log(p_{j|i})`} ` 34:T8a6,{ const dist_and_probs = anchor_distances .filter(d => d.distance !== 0) .map(d => ({dist: d.distance, prob:gaussian_prob(d.distance, variance)})), distance_max = d3.max(anchor_distances, d => d.distance), entropy = calc_entropy(dist_and_probs.map(d => d.prob)), prob_dist_curve = linspace(0, distance_max, 300) .map(dist => ({dist, prob:gaussian_prob(dist, variance)})), prob_extents = d3.extent(prob_dist_curve, d => d.prob); const x = d3.scaleLinear() .domain([0, distance_max]) .range([margin.left, width - margin.right]), y = d3.scaleLinear() .domain(prob_extents) .range([height - margin.bottom, margin.top]), point_color = d3.scaleLinear() .domain(prob_extents) .range(['grey','steelblue']), area = d3.area() .x(d => x(d.dist)) .y0(y(0)) .y1(d => y(d.prob)); const svg = d3.select(DOM.svg(width, height)); svg.append("text") .text(`Entropy: ${entropy.toPrecision(3)}`) .attr('x', width - 10) .attr('y', 40) .attr('text-anchor', 'end') .attr('font-size', 20) svg.append("text") .text(`Perplexity: ${exp(entropy).toPrecision(3)}`) .attr('x', width - 10) .attr('y', 70) .attr('text-anchor', 'end') .attr('font-size', 20) svg.append("text") .text(`Perplexity Goal: ${perplexity.toPrecision(3)}`) .attr('x', width - 10) .attr('y', 100) .attr('text-anchor', 'end') .attr('font-size', 20) svg.append("path") .datum(prob_dist_curve) .attr("fill", "steelblue") .attr("fill-opacity", 0.2) .attr("d", area); svg.append("circle") .attr("cx", x(0)) .attr("cy", y(0)) .attr("stroke", "grey") .attr("stroke-opacity", 0.2) .attr('fill', 'orangered') .attr('fill-opacity', 0.75) .attr("r", 7); svg.append("g") .selectAll("line") .data(dist_and_probs) .enter().append("line") .attr("x1", d => x(d.dist)) .attr("x2", d => x(d.dist)) .attr("y1", y(0) - 10) .attr("y2", y(0) + 10) .attr("stroke", d => point_color(d.prob)) .attr("stroke-width", 2) .attr("stroke-opacity", 0.5); return svg.node(); }35:T618,function variance_binary_search(i, pw_dists, entropy_tol = 1e-3, max_attempts = 40){ // grab the distances for point i const dists_from_i = pw_dists .reduce( (dists, d) => d.source === i ? [...dists, d.distance] : dists, [] ); let variance_min = 0, // lower bound variance_max = Infinity, // upper bound variance = 1, // initial value searching = true, attempts = 0, probs_given_variance, current_entropy; // number of iterations through while looking while(searching){ attempts++ // increment out attempts counter // calculate the probabilities and entropy given the current beta. probs_given_variance = calc_gaussian_probs(dists_from_i, variance); current_entropy = calc_entropy(probs_given_variance); // Entropy is below where we want, so try a higher variance if(current_entropy < target_entropy){ variance_min = variance; // move our min bound up variance = variance_max === Infinity ? variance*2 : (variance+variance_max)/2; } else { // distribution is too narrow. variance_max = variance; // move our max bound down variance = variance_min === -Infinity ? variance/2 : (variance+variance_min)/2; } // check some termination conditions const too_many_attempts = attempts > max_attempts; const close_to_target = Math.abs(current_entropy - target_entropy) < entropy_tol; searching = !(too_many_attempts || close_to_target); } return probs_given_variance }36:T4ab,md` # The mapping calculations We have now done all we need to do with the high dimensional data. Since at each iteration the high dimensional probabilities don't move we only have to do everything above once. From now on, what we need to do is take a random initial mapping of all our datapoints into two(or three) dimensions: ${tex`y_{\{1,2, ..., ${num_points}\}}`}, and then nudge those mappings around until they best 'match' the distribution of probabilities we generated in the above sections. What 'match' means in this context we will get into in two seconds, but first we need to be able to take our mapped distances and translate them to probabilities. Just like we did in the high dimension. ## Calculating the mapping probabilities We will first make a function that calculates the conditional probabilities in the low-dimensional mapping. This function looks very similar to the one for the high dimensional mapings, but now it is the [t-distribution](https://en.wikipedia.org/wiki/Student%27s_t-distribution) with one degree of freedom (aka the [Cauchy](https://en.wikipedia.org/wiki/Cauchy_distribution) distribution). We will denote these probabilties with ${tex`q_{ij}`}. `37:T46f,md`Note that we don't have any ${tex`\sigma`} to tune like we did in the high-dimension case, as we simply set it to a constant and it cancels out. This also contributes to the fact that our normalization constant (the denomonator) is over all pairs of points, not just for a common point. Another important difference is there is no exponential in the formula. This is nice because exponentials are expensive to calculate. ## The cost function Next, we need some form of telling the model how good it is doing. What we want is our probabilities from the high dimension ${tex`P`}to match the probabilities for the low dimension ${tex`Q`} as much as possible. Luckily in statistics there is an equation for doing just this. ### Kullback-Leibler Divergence The [Kullback-Leibler Divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence) (or KLD) is a measure of how close two distributions are to eachother. It's not a distance in the strict definition but instead a similarity measure. It's used all over in machine learning and statistics to measure goodness of fit. The equation is fortunately very simple.`38:T428,function step_down_gradient(iter, Y, Y_step){ const momentum = momentums[iter > momentum_switch ? 1: 0], Q = find_Q(Y); let average_y = make_array(map_dims).map(d=>0), cost_total = 0, grad_id, step_id,step_id_new; for(let i=0; i y[0]); const dim2_domain = d3.extent(Y, (y) => y[1]); const dim1 = d3.scaleLinear().domain(dim1_domain).range([margin.left, width - margin.right]); const dim2 = d3.scaleLinear().domain(dim2_domain).range([margin.top, height - margin.bottom]); const group_color = d3.scaleOrdinal() .domain(make_array(num_clusters).map((d,i)=>i)) .range(d3_color.schemeCategory10); const vis_data = data.map((d,i) => ({...d, map: Y[i]})); svg.append('text') .text(`Iteration: ${costs.length}`) .attr('x', margin.right) .attr('y', 2*margin.top) .attr('font-size', 30) .attr('font-family', 'optima') svg.append("g") .attr("stroke", "#000") .attr("stroke-opacity", 0.2) .selectAll("circle") .data(vis_data) .enter().append("circle") .attr("cx", d => dim1(d.map[0])) .attr("cy", d => dim2(d.map[1])) .attr("fill", d => group_color(d.group)) .attr("r", 5); return svg.node(); }3b:T646,{ const x = d3.scaleLinear() .domain([0, costs.length]) .range([margin.left, width - margin.right]), y = d3.scaleLinear() .domain(d3.extent(costs)) .range([height - margin.bottom, margin.top]), line = d3.line() .x((d,i) => x(i)) .y((d) => y(d)); const svg = d3.select(DOM.svg(width, height)); svg.append('line') .attr('x1', x(early_exag_itts)) .attr('x2', x(early_exag_itts)) .attr('y1', height - margin.top) .attr('y2', margin.bottom) .attr('stroke', 'steelblue') .attr('stroke-width', 2) svg.append('text') .text('Early exageration switch off') .attr('x', x(early_exag_itts) - 5) .attr('y', 3*margin.top) .attr('font-family', 'optima') .attr('text-anchor', 'end') svg.append('line') .attr('x1', x(momentum_switch)) .attr('x2', x(momentum_switch)) .attr('y1', height - margin.top) .attr('y2', margin.bottom) .attr('stroke', 'steelblue') .attr('stroke-width', 2) svg.append('text') .text(`Momentum changed to ${momentum_2}`) .attr('x', x(momentum_switch) + 5) .attr('y', 5*margin.top) .attr('font-family', 'optima') svg.append('text') .attr('x', width - 50) .attr('y', 50) .attr('text-anchor', 'end') .attr('font-size', 35) .attr('font-family', 'optima') .text('Cost'); svg.append("path") .datum(costs) .attr("fill", "none") .attr("stroke", "orangered") .attr("stroke-width", 2) .attr("stroke-linejoin", "round") .attr("stroke-linecap", "round") .attr("d", line); return svg.node(); }3c:T424,data = { const make_centroid = (ndim) => { 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) }3d:[] 2d:{"title":"t-SNE explained in plain javascript","cells":[{"id":0,"value":"md`# t-SNE explained in plain javascript`","pinned":false,"hidden":false,"mode":"ojs"},{"id":54,"value":"$2f","pinned":false,"hidden":false,"mode":"ojs"},{"id":835,"value":"viewof n = slider3({\n message:'Set number of observations in your dataset',\n min: 1, \n max: 100,\n start: 50,\n step: 1, \n int:true})","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":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":1984,"value":"calc_distance = (x1, x2) => sqrt(x1.reduce(\n (sum, x1_i, i) => sum + squared(x1_i - x2[i]),\n 0\n))","pinned":false,"hidden":false,"mode":"ojs"},{"id":895,"value":"//Get index in a 1-d array from a 2d matrix reference point. (Assumes our matrix is square)\n//This function will auto update when we change the number of points so we don't need to pass \n//that as a parameter thanks to the beauty of reactive programing!\nmat_ind = (i,j) => i*num_points + j ","pinned":false,"hidden":false,"mode":"ojs"},{"id":1987,"value":"squared = (val) => Math.pow(val, 2)","pinned":false,"hidden":false,"mode":"ojs"},{"id":906,"value":"md`Okay, now we can get into the pairwise distances calculation itself. `","pinned":false,"hidden":false,"mode":"ojs"},{"id":62,"value":"pairwise_distances = {\n // initialize some variables we use inside our loops\n let x_i, x_j, distance;\n // make_array is another helper function that just generates a new empty array of a given size.\n const pw_dists = make_array(squared(num_points));\n for(let i = 0; i < num_points; i++){\n \tx_i = data[i].location;\n // deal with the point itself here.\n pw_dists[mat_ind(i,i)] = {source: i, target: i, distance:0};\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[mat_ind(i,j)] = {source: i, target: j, distance};\n pw_dists[mat_ind(j,i)] = {source: j, target: i, distance};\n }\n }\n return pw_dists\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":59,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":918,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":967,"value":"viewof new_point = slider3(\n {message:'Set distance of point j from point i',\n min:-5, max:5, \n start:1.1, step: 0.1, precision: 2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1022,"value":"viewof orig_variance = slider3(\n {message:'Set the variance of the distribution',\n min:0.1, max:3, \n start:1, step: 0.1, precision: 2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1062,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":2033,"value":"function calc_entropy(probs){\n return probs.reduce( \n (sum, p) => sum - (p>1e-7 ? p*log(p): 0), // super small probabilities can be ignored\n 0\n )\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2032,"value":"\nmd`\nBelow we plot a point in our data and each other point's distance away from that point. The gaussian probability curve is drawn over these points based upon the current variance for the distribution drawn above the given point. As you change the variance with the slider notice how the distribution gets wider or smaller to essentially 'grab' more or less points that contribute a meaningful amount of probability to it's overall sum of probabilities. \n\n\n`\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":1250,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":1503,"value":"viewof variance = slider3({\n message: 'Set the variance for above', \n min: 1,start:2, max: 50, step: 0.1, precision: 3})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1268,"value":"anchor_distances = {\n const anchor_point = 5;\n const anchor_distances = pairwise_distances\n \t.filter(d => d.source === anchor_point);\n return anchor_distances\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1255,"value":"md`\nIn order to find the optimal variance for each point in our dataset we use a searching technique called [binary search](https://en.wikipedia.org/wiki/Binary_search_algorithm). What binary search is doing here is goes point by point and searching in a disciplined way for the variance value that gets the perplexity as close to the desired perplexity as possible. Just like you did with the slider above.\n\nTo do this we first need to find our goal entropy for our points. Do do this we simply take the log of our goal perplexity (in this case __${perplexity}__).\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1274,"value":"target_entropy = log(perplexity) // entropy level we're shooting for given perplexity","pinned":true,"hidden":false,"mode":"ojs"},{"id":1278,"value":"md`\nNow that we have the entropy we just need a couple of functions to make our lives easier.\n\nThe first is one that will calculate the un-normalized gaussian probability of point \\`j\\` given its distance from point \\`i\\` and point \\`i\\`'s variance...\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2371,"value":"tex.block`p_{j|i} \\propto \\exp\\{-(\\text{distance from } i \\text{ to } j)^2 / 2\\sigma_i\\}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1095,"value":"function gaussian_prob(distance, variance){\n return exp(-squared(distance) / (2*variance))\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1292,"value":"md`\nNext is a function that takes an array of distances from point \\`i\\` and converts them to an array containing of normalized probabilities according to the guassian equation above.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2367,"value":"tex.block`\nP_{i|j} = \\frac{\\exp(-[\\text{distance between } i \\text{ and } j]^2 / 2\\sigma_i)}{\\sum_{k \\ne i} \\exp(-[\\text{distance between } i \\text{ and } k]^2 / 2\\sigma_i)}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1429,"value":"md`Essentially for every combination we have the un-normalized probability of \\`j\\` given \\`i\\` as the anchor, which is then divided by the sum of all the probabilities given point \\`i\\` to normalize it.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1170,"value":"function calc_gaussian_probs(distances, variance){\n let sum_probs = 0;\n\n // calc un-normalized probs and accumulate sum of probs for normalization\n return distances.map( (dist, i) => {\n const curr_prob = dist !== 0 ? gaussian_prob(dist, variance): 0;\n sum_probs += curr_prob; // accumulate the sums \t\n return curr_prob; // send probability to our array \n }).map(p => p/sum_probs) // return the prob vector after normalizing it using sum.\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1071,"value":"md`\nLastly, we set up a function for doing the binary search for the optimal variance for a given point. I won't go too far into this other than to leave the code here. What matters is that the algorithm bounces around checking the perplexity that results from different variance values in a disciplined way in an effort to make the entropy match our goal, it then returns the vector of probabilities that correspond to the optimal variance.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1531,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":273,"value":"md`\n## Calculating probabilities\n\nNow that we have our binary search function in place, we can calculate the non-symmetric probabilities. We do this for each datapoint to build out our pairwise probabilties for every pairwise connection between the points. \n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":108,"value":"P_i_given_j = {\n let p_row, i, j;\n const P_i_given_j = make_array(num_points*num_points); \n \n for(i = 0; i < num_points; i++){\n p_row = variance_binary_search(i, pairwise_distances) // probs from best variance. \n for(j=0;j q/q_sum);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1614,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":2381,"value":"tex.block`C = \\text{KLD}(P,Q) = \\sum_i\\sum_jp_{ij}\\log\\left(\\frac{p_{ij}}{q_{ij}}\\right)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2385,"value":"md`Below we will define a function for performing the inside of those two sums which we will then use to accumulate the full KLD when running our algorithm.\n\n_In an earlier version of this I miswrote the KLD measure here: putting a negative sign in front of it. This was kindly pointed out [on twitter by Mark Probst](https://twitter.com/markprobst). I had mixed up the definition of the cost function and the goal of minimizing the cost function. In addition the actual code for it was wrong and I forgot to put the fraction inside of the \\`log\\` and instead left just \\`q\\`. Oops!_`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1693,"value":"function KLD(p, q){\n return p * log(p/q)\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1691,"value":"md`\n## The gradient\nNow that we have our cost function set up the one thing remaining is a way to calculate the gradient of our mapping positions ${tex`\\underline{Y}`} with respect to that cost function. Like the KLD formula, we are lucky here in that it is an very simple (relatively) formula:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2391,"value":"tex.block`\\frac{\\partial C}{\\partial y_i} = 4\\sum_j\\left[(p_{ij} - q_{ij})(y_i - y_j)(1 + (\\text{distance from }y_i \\text{ to }y_j)^2)^{-1}\\right]`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2393,"value":"md`I won't bore you with the math behind it but just know that this says, \"given we have a point's mapped position (${tex`y_i`}) along with its corresponding high dimensional probabilities (${tex`p_{i\\cdot}`}) and low dimension probabilities (${tex`q_{i\\cdot}`}), we can figure out what direction to nudge ${tex`y_i`} in order to lower the cost ${tex`C`} of our mapping.\"\n\nWe take advantage of the fact that the same loop can be used to accumulate both the cost and the gradient and define a function that returns both for a given point.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1745,"value":"function cost_grad_calc(Y, i, P, Q, iter){\n let grad = make_array(map_dims).map( () => 0),\n cost = 0,\n // the paper recomends exagerating the distances of the high dimensional points\n // for a few iterations at first, this checks if we are still exagerating.\n early_exag = iter < early_exag_itts? early_exag_amount: 1,\n grad_scaler, p_ij, q_ij;\n \n for(let j=0; j`\n \n div.style = \"border-radius:5px; width: 100px; height: 30px; color: white; background-color: orangered;\";\n return div\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":581,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs"},{"id":729,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs"},{"id":2132,"value":"md`__Gradient descent parameters__`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2323,"value":"viewof perplexity = slider3({\n message: \"Perplexity value\",\n min: 1, max: 60, start: 12,\n step: 1, int: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":770,"value":"viewof learning_rate = slider3({\n message: \"Learning rate\",\n min: 0.001, max: 20, start: 2,\n step: 0.001})","pinned":false,"hidden":false,"mode":"ojs"},{"id":725,"value":"viewof num_itterations = slider3({\n message: \"Number of iterations to run algorithm\",\n min: 1, max: 500, start: 200,\n step: 1, int: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2134,"value":"viewof early_exag_itts = slider3({\n message: \"How many iterations to exagerate\",\n min: 0, max: num_itterations, start: 75,\n step: 1, int: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":574,"value":"viewof early_exag_amount = slider3({\n message: \"How much to exagerate high dimensional distances at first\",\n min: 1, max: 10, start: 4,\n step: 0.5})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1847,"value":"viewof momentum_switch = slider3({\n message: \"How many iterations until switching over momentum amount\",\n min: 0, max: num_itterations, start: Math.round(num_itterations/2),\n step: 1, int: true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1855,"value":"viewof momentum_1 = slider3({\n message: \"First momentum amount\",\n start: 0.5, step : 0.05, precision: 2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1859,"value":"viewof momentum_2 = slider3({\n message: \"Second momentum amount\",\n start: 0.8, step : 0.05, precision: 2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2095,"value":"md`\n__Data Generation__\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2097,"value":"viewof num_dimensions = slider3({\n message: \"Number of dimensions in original data\",\n start:10, min:2, max:15, step : 1, precision: 1, int:true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2102,"value":"viewof num_clusters = slider3({\n message: \"Number of clusters in data\",\n start: 5, min:1, max:10, step : 1, precision: 1, int:true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2106,"value":"viewof num_samples = slider3({\n message: \"Number of samples drawn from each cluster\",\n start: 25, min:1, max:100, step : 1, precision: 1, int:true})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2116,"value":"viewof cluster_spread = slider3({\n message: \"How spread apart clusters can be in each dimension\",\n start: 15, min:0, max:20, step : 0.5, precision: 2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2123,"value":"viewof random_noise_sd = slider3({\n message: \"How dispursed each cluster is\",\n start: 2, min:0, max:10, step : 0.5, precision: 2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":2187,"value":"md`__Helper functions and various sundries__\n\nBelow are an assortment of intermediate variables that are automatically calculated and helper functions I have written (mostly just removing the \\`Math.\\` from basic math operations.)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":483,"value":"map_dims = 2","pinned":false,"hidden":false,"mode":"ojs"},{"id":2145,"value":"momentums = [momentum_1, momentum_2]","pinned":false,"hidden":false,"mode":"ojs"},{"id":864,"value":"num_points = data.length","pinned":false,"hidden":false,"mode":"ojs"},{"id":140,"value":"exp = Math.exp","pinned":false,"hidden":false,"mode":"ojs"},{"id":1458,"value":"sqrt = Math.sqrt","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":1910,"value":"function initialize_points({n, dim = 2, fill}){\n const fill_function = typeof fill == 'undefined' ?\n Math.random :\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":932,"value":"function linspace(low, high, n = 100){\n const width = (high - low), half = width/2;\n return make_array(100).map((d,i) => low + ((i/100)*width))\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2213,"value":"function slider3({\n message = 'Set your value', \n min = 0, max = 1, start = 0.5, step = 0.1,\n precision = 3, int = false,\n}){\n \n min = min/step;\n max = max/step;\n start = start/step;\n \n const div = html`\n\t${message}\n \n ${start*step}\n `;\n const range = div.querySelector(\"[type=range]\");\n const value_disp = div.querySelector('span');\n div.value = value_disp.innerHTML = start*step;\n\n \n range.addEventListener(\"input\", () => {\n const newValue = range.valueAsNumber*step; \n div.value = newValue;\n value_disp.innerHTML = int ? Math.round(newValue) : newValue.toPrecision(precision); \n });\n return div;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2229,"value":"md`This wonderful function for combining markdown and tex into the same cell was taken from [this notebook.](https://beta.observablehq.com/@kelleyvanevert/katex-within-markdown) `","pinned":false,"hidden":false,"mode":"ojs"},{"id":2200,"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":2203,"value":"$3c","pinned":false,"hidden":false,"mode":"ojs"},{"id":2210,"value":"md`Variables for charts`","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":2224,"value":"md`__Libraries!__`","pinned":false,"hidden":false,"mode":"ojs"},{"id":721,"value":"d3_color = require('d3-scale-chromatic')","pinned":false,"hidden":false,"mode":"ojs"},{"id":47,"value":"d3 = require('d3')","pinned":false,"hidden":false,"mode":"ojs"},{"id":915,"value":"stats = require('statdists')","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W3d","stdlib":"1","titling":"auto"}