qdtsne
Quick and dirty t-SNE in C++
Loading...
Searching...
No Matches
Status.hpp
Go to the documentation of this file.
1/*
2 *
3 * Copyright (c) 2014, Laurens van der Maaten (Delft University of Technology)
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions are met:
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 * 3. All advertising materials mentioning features or use of this software
14 * must display the following acknowledgement:
15 * This product includes software developed by the Delft University of Technology.
16 * 4. Neither the name of the Delft University of Technology nor the names of
17 * its contributors may be used to endorse or promote products derived from
18 * this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY LAURENS VAN DER MAATEN ''AS IS'' AND ANY EXPRESS
21 * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
22 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
23 * EVENT SHALL LAURENS VAN DER MAATEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
24 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
26 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
28 * IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
29 * OF SUCH DAMAGE.
30 *
31 */
32
33#ifndef QDTSNE_STATUS_HPP
34#define QDTSNE_STATUS_HPP
35
36#include <vector>
37#include <array>
38#include <algorithm>
39#include <cstddef>
40
41#include "sanisizer/sanisizer.hpp"
42
43#include "SPTree.hpp"
44#include "Options.hpp"
45#include "utils.hpp"
46
52namespace qdtsne {
53
64template<std::size_t num_dim_, typename Index_, typename Float_>
65class Status {
66public:
70 Status(NeighborList<Index_, Float_> neighbors, Options options) :
71 my_neighbors(std::move(neighbors)),
72 my_tree(sanisizer::cast<internal::SPTreeIndex>(my_neighbors.size()), options.max_depth),
73 my_options(std::move(options))
74 {
75 const auto nobs = sanisizer::cast<Index_>(my_neighbors.size()); // use Index_ to check safety of cast in num_observations().
76 const std::size_t buffer_size = sanisizer::product<std::size_t>(nobs, num_dim_);
77
78 sanisizer::resize(my_dY, buffer_size);
79 sanisizer::resize(my_uY, buffer_size);
80 sanisizer::resize(my_gains, buffer_size, static_cast<Float_>(1));
81 sanisizer::resize(my_pos_f, buffer_size);
82 sanisizer::resize(my_neg_f, buffer_size);
83
84 if (options.num_threads > 1) {
85 sanisizer::resize(my_parallel_buffer, nobs);
86 }
87 }
92private:
93 NeighborList<Index_, Float_> my_neighbors;
94 std::vector<Float_> my_dY, my_uY, my_gains, my_pos_f, my_neg_f;
95
96 internal::SPTree<num_dim_, Float_> my_tree;
97 std::vector<Float_> my_parallel_buffer; // Buffer to hold parallel-computed results prior to reduction.
98
99 Options my_options;
100 int my_iter = 0;
101
102 typename decltype(I(my_tree))::LeafApproxWorkspace my_leaf_workspace;
103
104public:
108 int iteration() const {
109 return my_iter;
110 }
111
115 int max_iterations() const {
116 return my_options.max_iterations;
117 }
118
122 Index_ num_observations() const {
123 return my_neighbors.size(); // safety of the cast to Index_ is already checked in the constructor.
124 }
125
126#ifndef NDEBUG
130 const auto& get_neighbors() const {
131 return my_neighbors;
132 }
136#endif
137
138public:
150 void run(Float_* const Y, const int limit) {
151 Float_ multiplier = (my_iter < my_options.early_exaggeration_iterations ? my_options.exaggeration_factor : 1);
152 Float_ momentum = (my_iter < my_options.momentum_switch_iterations ? my_options.start_momentum : my_options.final_momentum);
153
154 for(; my_iter < limit; ++my_iter) {
155 if (my_iter == my_options.early_exaggeration_iterations) {
156 multiplier = 1;
157 }
158 if (my_iter == my_options.momentum_switch_iterations) {
159 momentum = my_options.final_momentum;
160 }
161 iterate(Y, multiplier, momentum);
162 }
163 }
164
174 void run(Float_* const Y) {
175 run(Y, my_options.max_iterations);
176 }
177
178private:
179 static Float_ sign(const Float_ x) {
180 constexpr Float_ zero = 0;
181 constexpr Float_ one = 1;
182 return (x == zero ? zero : (x < zero ? -one : one));
183 }
184
185 void iterate(Float_* const Y, const Float_ multiplier, const Float_ momentum) {
186 compute_gradient(Y, multiplier);
187
188 // Update gains
189 const auto buffer_size = my_gains.size();
190 for (decltype(I(buffer_size)) i = 0; i < buffer_size; ++i) {
191 Float_& g = my_gains[i];
192 constexpr Float_ lower_bound = 0.01;
193 constexpr Float_ to_add = 0.2;
194 constexpr Float_ to_mult = 0.8;
195 g = std::max(lower_bound, sign(my_dY[i]) != sign(my_uY[i]) ? (g + to_add) : (g * to_mult));
196 }
197
198 // Perform gradient update (with momentum and gains)
199 for (decltype(I(buffer_size)) i = 0; i < buffer_size; ++i) {
200 my_uY[i] = momentum * my_uY[i] - my_options.eta * my_gains[i] * my_dY[i];
201 Y[i] += my_uY[i];
202 }
203
204 // Make solution zero-mean for each dimension
205 std::array<Float_, num_dim_> means{};
206 const Index_ num_obs = num_observations();
207 for (Index_ i = 0; i < num_obs; ++i) {
208 for (std::size_t d = 0; d < num_dim_; ++d) {
209 means[d] += Y[sanisizer::nd_offset<std::size_t>(d, num_dim_, i)];
210 }
211 }
212 for (std::size_t d = 0; d < num_dim_; ++d) {
213 means[d] /= num_obs;
214 }
215 for (Index_ i = 0; i < num_obs; ++i) {
216 for (std::size_t d = 0; d < num_dim_; ++d) {
217 Y[sanisizer::nd_offset<std::size_t>(d, num_dim_, i)] -= means[d];
218 }
219 }
220
221 return;
222 }
223
224private:
225 void compute_gradient(const Float_* const Y, const Float_ multiplier) {
226 my_tree.set(Y);
227 compute_edge_forces(Y, multiplier);
228
229 std::fill(my_neg_f.begin(), my_neg_f.end(), 0);
230 const Float_ sum_Q = compute_non_edge_forces();
231
232 const auto buffer_size = my_dY.size();
233 for (decltype(I(buffer_size)) i = 0; i < buffer_size; ++i) {
234 my_dY[i] = my_pos_f[i] - (my_neg_f[i] / sum_Q);
235 }
236 }
237
238 void compute_edge_forces(const Float_* const Y, Float_ multiplier) {
239 std::fill(my_pos_f.begin(), my_pos_f.end(), 0);
240 const Index_ num_obs = num_observations();
241
242 parallelize(my_options.num_threads, num_obs, [&](const int, const Index_ start, const Index_ length) -> void {
243 for (Index_ i = start, end = start + length; i < end; ++i) {
244 const auto& current = my_neighbors[i];
245 const auto offset = sanisizer::product_unsafe<std::size_t>(i, num_dim_);
246 const auto self = Y + offset;
247 const auto pos_out = my_pos_f.data() + offset;
248
249 for (const auto& x : current) {
250 Float_ sqdist = 0;
251 const auto neighbor = Y + sanisizer::product_unsafe<std::size_t>(x.first, num_dim_);
252 for (std::size_t d = 0; d < num_dim_; ++d) {
253 Float_ delta = self[d] - neighbor[d];
254 sqdist += delta * delta;
255 }
256
257 const Float_ mult = multiplier * x.second / (static_cast<Float_>(1) + sqdist);
258 for (std::size_t d = 0; d < num_dim_; ++d) {
259 pos_out[d] += mult * (self[d] - neighbor[d]);
260 }
261 }
262 }
263 });
264
265 return;
266 }
267
268 Float_ compute_non_edge_forces() {
269 if (my_options.leaf_approximation) {
270 my_tree.compute_non_edge_forces_for_leaves(my_options.theta, my_leaf_workspace, my_options.num_threads);
271 }
272
273 const Index_ num_obs = num_observations();
274 if (my_options.num_threads > 1) {
275 // Don't use reduction methods, otherwise we get numeric imprecision
276 // issues (and stochastic results) based on the order of summation.
277 parallelize(my_options.num_threads, num_obs, [&](const int, const Index_ start, const Index_ length) -> void {
278 for (Index_ i = start, end = start + length; i < end; ++i) {
279 const auto neg_ptr = my_neg_f.data() + sanisizer::product_unsafe<std::size_t>(i, num_dim_);
280 if (my_options.leaf_approximation) {
281 my_parallel_buffer[i] = my_tree.compute_non_edge_forces_from_leaves(i, neg_ptr, my_leaf_workspace);
282 } else {
283 my_parallel_buffer[i] = my_tree.compute_non_edge_forces(i, my_options.theta, neg_ptr);
284 }
285 }
286 });
287 return std::accumulate(my_parallel_buffer.begin(), my_parallel_buffer.end(), static_cast<Float_>(0));
288 }
289
290 Float_ sum_Q = 0;
291 for (Index_ i = 0; i < num_obs; ++i) {
292 const auto neg_ptr = my_neg_f.data() + sanisizer::product_unsafe<std::size_t>(i, num_dim_);
293 if (my_options.leaf_approximation) {
294 sum_Q += my_tree.compute_non_edge_forces_from_leaves(i, neg_ptr, my_leaf_workspace);
295 } else {
296 sum_Q += my_tree.compute_non_edge_forces(i, my_options.theta, neg_ptr);
297 }
298 }
299 return sum_Q;
300 }
301};
302
303}
304
305#endif
Options for the t-SNE algorithm.
Status of the t-SNE iterations.
Definition Status.hpp:65
void run(Float_ *const Y, const int limit)
Definition Status.hpp:150
int iteration() const
Definition Status.hpp:108
Index_ num_observations() const
Definition Status.hpp:122
void run(Float_ *const Y)
Definition Status.hpp:174
int max_iterations() const
Definition Status.hpp:115
Quick and dirty t-SNE.
knncolle::NeighborList< Index_, Float_ > NeighborList
Lists of neighbors for each observation.
Definition utils.hpp:37
void parallelize(const int num_workers, const Task_ num_tasks, Run_ run_task_range)
Definition utils.hpp:124
Options for initialize().
Definition Options.hpp:14
int max_depth
Definition Options.hpp:105
int momentum_switch_iterations
Definition Options.hpp:73
int early_exaggeration_iterations
Definition Options.hpp:58
int max_iterations
Definition Options.hpp:50
bool leaf_approximation
Definition Options.hpp:117
double eta
Definition Options.hpp:91
double final_momentum
Definition Options.hpp:85
double exaggeration_factor
Definition Options.hpp:64
int num_threads
Definition Options.hpp:124
double theta
Definition Options.hpp:44
Utilities for running t-SNE.