Faiss

- Main Page
- Namespaces
- Classes
- Files

- File List

All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends
/data/users/matthijs/github_faiss/faiss/VectorTransform.cpp
1 /**
2 * Copyright (c) 2015-present, Facebook, Inc.
3 * All rights reserved.
4 *
5 * This source code is licensed under the BSD+Patents license found in the
6 * LICENSE file in the root directory of this source tree.
7 */
8
9 // Copyright 2004-present Facebook. All Rights Reserved
10 // -*- c++ -*-
11
12 #include "VectorTransform.h"
13
14 #include <cstdio>
15 #include <cmath>
16 #include <cstring>
17
18 #include "utils.h"
19 #include "FaissAssert.h"
20 #include "IndexPQ.h"
21
22 using namespace faiss;
23
24
25 extern "C" {
26
27 // this is to keep the clang syntax checker happy
28 #ifndef FINTEGER
29 #define FINTEGER int
30 #endif
31
32
33 /* declare BLAS functions, see http://www.netlib.org/clapack/cblas/ */
34
35 int sgemm_ (
36 const char *transa, const char *transb, FINTEGER *m, FINTEGER *
37 n, FINTEGER *k, const float *alpha, const float *a,
38 FINTEGER *lda, const float *b,
39 FINTEGER *ldb, float *beta,
40 float *c, FINTEGER *ldc);
41
42 int ssyrk_ (
43 const char *uplo, const char *trans, FINTEGER *n, FINTEGER *k,
44 float *alpha, float *a, FINTEGER *lda,
45 float *beta, float *c, FINTEGER *ldc);
46
47 /* Lapack functions from http://www.netlib.org/clapack/old/single/ */
48
49 int ssyev_ (
50 const char *jobz, const char *uplo, FINTEGER *n, float *a,
51 FINTEGER *lda, float *w, float *work, FINTEGER *lwork,
52 FINTEGER *info);
53
54 int dsyev_ (
55 const char *jobz, const char *uplo, FINTEGER *n, double *a,
56 FINTEGER *lda, double *w, double *work, FINTEGER *lwork,
57 FINTEGER *info);
58
59 int sgesvd_(
60 const char *jobu, const char *jobvt, FINTEGER *m, FINTEGER *n,
61 float *a, FINTEGER *lda, float *s, float *u, FINTEGER *ldu, float *vt,
62 FINTEGER *ldvt, float *work, FINTEGER *lwork, FINTEGER *info);
63
64 }
65
66 /*********************************************
67 * VectorTransform
68 *********************************************/
69
70
71
72 float * VectorTransform::apply (Index::idx_t n, const float * x) const
73 {
74 float * xt = new float[n * d_out];
75 apply_noalloc (n, x, xt);
76 return xt;
77 }
78
79
80 void VectorTransform::train (idx_t, const float *) {
81 // does nothing by default
82 }
83
84
85 void VectorTransform::reverse_transform (
86 idx_t , const float *,
87 float *) const
88 {
89 FAISS_THROW_MSG ("reverse transform not implemented");
90 }
91
92
93
94
95 /*********************************************
96 * LinearTransform
97 *********************************************/
98
99 /// both d_in > d_out and d_out < d_in are supported
100 LinearTransform::LinearTransform (int d_in, int d_out,
101 bool have_bias):
102 VectorTransform (d_in, d_out), have_bias (have_bias),
103 is_orthonormal (false), verbose (false)
104 {}
105
106 void LinearTransform::apply_noalloc (Index::idx_t n, const float * x,
107 float * xt) const
108 {
109 FAISS_THROW_IF_NOT_MSG(is_trained, "Transformation not trained yet");
110
111 float c_factor;
112 if (have_bias) {
113 FAISS_THROW_IF_NOT_MSG (b.size() == d_out, "Bias not initialized");
114 float * xi = xt;
115 for (int i = 0; i < n; i++)
116 for(int j = 0; j < d_out; j++)
117 *xi++ = b[j];
118 c_factor = 1.0;
119 } else {
120 c_factor = 0.0;
121 }
122
123 FAISS_THROW_IF_NOT_MSG (A.size() == d_out * d_in,
124 "Transformation matrix not initialized");
125
126 float one = 1;
127 FINTEGER nbiti = d_out, ni = n, di = d_in;
128 sgemm_ ("Transposed", "Not transposed",
129 &nbiti, &ni, &di,
130 &one, A.data(), &di, x, &di, &c_factor, xt, &nbiti);
131
132 }
133
134
135 void LinearTransform::transform_transpose (idx_t n, const float * y,
136 float *x) const
137 {
138 if (have_bias) { // allocate buffer to store bias-corrected data
139 float *y_new = new float [n * d_out];
140 const float *yr = y;
141 float *yw = y_new;
142 for (idx_t i = 0; i < n; i++) {
143 for (int j = 0; j < d_out; j++) {
144 *yw++ = *yr++ - b [j];
145 }
146 }
147 y = y_new;
148 }
149
150 {
151 FINTEGER dii = d_in, doi = d_out, ni = n;
152 float one = 1.0, zero = 0.0;
153 sgemm_ ("Not", "Not", &dii, &ni, &doi,
154 &one, A.data (), &dii, y, &doi, &zero, x, &dii);
155 }
156
157 if (have_bias) delete [] y;
158 }
159
160 void LinearTransform::set_is_orthonormal ()
161 {
162 if (d_out > d_in) {
163 // not clear what we should do in this case
164 is_orthonormal = false;
165 return;
166 }
167 if (d_out == 0) { // borderline case, unnormalized matrix
168 is_orthonormal = true;
169 return;
170 }
171
172 double eps = 4e-5;
173 FAISS_ASSERT(A.size() >= d_out * d_in);
174 {
175 std::vector<float> ATA(d_out * d_out);
176 FINTEGER dii = d_in, doi = d_out;
177 float one = 1.0, zero = 0.0;
178
179 sgemm_ ("Transposed", "Not", &doi, &doi, &dii,
180 &one, A.data (), &dii,
181 A.data(), &dii,
182 &zero, ATA.data(), &doi);
183
184 is_orthonormal = true;
185 for (long i = 0; i < d_out; i++) {
186 for (long j = 0; j < d_out; j++) {
187 float v = ATA[i + j * d_out];
188 if (i == j) v-= 1;
189 if (fabs(v) > eps) {
190 is_orthonormal = false;
191 }
192 }
193 }
194 }
195
196 }
197
198
199 void LinearTransform::reverse_transform (idx_t n, const float * xt,
200 float *x) const
201 {
202 if (is_orthonormal) {
203 transform_transpose (n, xt, x);
204 } else {
205 FAISS_THROW_MSG ("reverse transform not implemented for non-orthonormal matrices");
206 }
207 }
208
209
210
211 /*********************************************
212 * RandomRotationMatrix
213 *********************************************/
214
215 void RandomRotationMatrix::init (int seed)
216 {
217
218 if(d_out <= d_in) {
219 A.resize (d_out * d_in);
220 float *q = A.data();
221 float_randn(q, d_out * d_in, seed);
222 matrix_qr(d_in, d_out, q);
223 } else {
224 A.resize (d_out * d_out);
225 float *q = A.data();
226 float_randn(q, d_out * d_out, seed);
227 matrix_qr(d_out, d_out, q);
228 // remove columns
229 int i, j;
230 for (i = 0; i < d_out; i++) {
231 for(j = 0; j < d_in; j++) {
232 q[i * d_in + j] = q[i * d_out + j];
233 }
234 }
235 A.resize(d_in * d_out);
236 }
237 is_orthonormal = true;
238 }
239
240 /*********************************************
241 * PCAMatrix
242 *********************************************/
243
244 PCAMatrix::PCAMatrix (int d_in, int d_out,
245 float eigen_power, bool random_rotation):
246 LinearTransform(d_in, d_out, true),
247 eigen_power(eigen_power), random_rotation(random_rotation)
248 {
249 is_trained = false;
250 max_points_per_d = 1000;
251 balanced_bins = 0;
252 }
253
254
255 namespace {
256
257 /// Compute the eigenvalue decomposition of symmetric matrix cov,
258 /// dimensions d_in-by-d_in. Output eigenvectors in cov.
259
260 void eig(size_t d_in, double *cov, double *eigenvalues, int verbose)
261 {
262 { // compute eigenvalues and vectors
263 FINTEGER info = 0, lwork = -1, di = d_in;
264 double workq;
265
266 dsyev_ ("Vectors as well", "Upper",
267 &di, cov, &di, eigenvalues, &workq, &lwork, &info);
268 lwork = FINTEGER(workq);
269 double *work = new double[lwork];
270
271 dsyev_ ("Vectors as well", "Upper",
272 &di, cov, &di, eigenvalues, work, &lwork, &info);
273
274 delete [] work;
275
276 if (info != 0) {
277 fprintf (stderr, "WARN ssyev info returns %d, "
278 "a very bad PCA matrix is learnt\n",
279 int(info));
280 // do not throw exception, as the matrix could still be useful
281 }
282
283
284 if(verbose && d_in <= 10) {
285 printf("info=%ld new eigvals=[", long(info));
286 for(int j = 0; j < d_in; j++) printf("%g ", eigenvalues[j]);
287 printf("]\n");
288
289 double *ci = cov;
290 printf("eigenvecs=\n");
291 for(int i = 0; i < d_in; i++) {
292 for(int j = 0; j < d_in; j++)
293 printf("%10.4g ", *ci++);
294 printf("\n");
295 }
296 }
297
298 }
299
300 // revert order of eigenvectors & values
301
302 for(int i = 0; i < d_in / 2; i++) {
303
304 std::swap(eigenvalues[i], eigenvalues[d_in - 1 - i]);
305 double *v1 = cov + i * d_in;
306 double *v2 = cov + (d_in - 1 - i) * d_in;
307 for(int j = 0; j < d_in; j++)
308 std::swap(v1[j], v2[j]);
309 }
310
311 }
312
313
314 }
315
316 void PCAMatrix::train (Index::idx_t n, const float *x)
317 {
318 const float * x_in = x;
319
320 x = fvecs_maybe_subsample (d_in, (size_t*)&n,
321 max_points_per_d * d_in, x, verbose);
322
323 ScopeDeleter<float> del_x (x != x_in ? x : nullptr);
324
325 // compute mean
326 mean.clear(); mean.resize(d_in, 0.0);
327 if (have_bias) { // we may want to skip the bias
328 const float *xi = x;
329 for (int i = 0; i < n; i++) {
330 for(int j = 0; j < d_in; j++)
331 mean[j] += *xi++;
332 }
333 for(int j = 0; j < d_in; j++)
334 mean[j] /= n;
335 }
336 if(verbose) {
337 printf("mean=[");
338 for(int j = 0; j < d_in; j++) printf("%g ", mean[j]);
339 printf("]\n");
340 }
341
342 if(n >= d_in) {
343 // compute covariance matrix, store it in PCA matrix
344 PCAMat.resize(d_in * d_in);
345 float * cov = PCAMat.data();
346 { // initialize with mean * mean^T term
347 float *ci = cov;
348 for(int i = 0; i < d_in; i++) {
349 for(int j = 0; j < d_in; j++)
350 *ci++ = - n * mean[i] * mean[j];
351 }
352 }
353 {
354 FINTEGER di = d_in, ni = n;
355 float one = 1.0;
356 ssyrk_ ("Up", "Non transposed",
357 &di, &ni, &one, (float*)x, &di, &one, cov, &di);
358
359 }
360 if(verbose && d_in <= 10) {
361 float *ci = cov;
362 printf("cov=\n");
363 for(int i = 0; i < d_in; i++) {
364 for(int j = 0; j < d_in; j++)
365 printf("%10g ", *ci++);
366 printf("\n");
367 }
368 }
369
370 std::vector<double> covd (d_in * d_in);
371 for (size_t i = 0; i < d_in * d_in; i++) covd [i] = cov [i];
372
373 std::vector<double> eigenvaluesd (d_in);
374
375 eig (d_in, covd.data (), eigenvaluesd.data (), verbose);
376
377 for (size_t i = 0; i < d_in * d_in; i++) PCAMat [i] = covd [i];
378 eigenvalues.resize (d_in);
379
380 for (size_t i = 0; i < d_in; i++)
381 eigenvalues [i] = eigenvaluesd [i];
382
383
384 } else {
385
386 std::vector<float> xc (n * d_in);
387
388 for (size_t i = 0; i < n; i++)
389 for(size_t j = 0; j < d_in; j++)
390 xc [i * d_in + j] = x [i * d_in + j] - mean[j];
391
392 // compute Gram matrix
393 std::vector<float> gram (n * n);
394 {
395 FINTEGER di = d_in, ni = n;
396 float one = 1.0, zero = 0.0;
397 ssyrk_ ("Up", "Transposed",
398 &ni, &di, &one, xc.data(), &di, &zero, gram.data(), &ni);
399 }
400
401 if(verbose && d_in <= 10) {
402 float *ci = gram.data();
403 printf("gram=\n");
404 for(int i = 0; i < n; i++) {
405 for(int j = 0; j < n; j++)
406 printf("%10g ", *ci++);
407 printf("\n");
408 }
409 }
410
411 std::vector<double> gramd (n * n);
412 for (size_t i = 0; i < n * n; i++)
413 gramd [i] = gram [i];
414
415 std::vector<double> eigenvaluesd (n);
416
417 // eig will fill in only the n first eigenvals
418
419 eig (n, gramd.data (), eigenvaluesd.data (), verbose);
420
421 PCAMat.resize(d_in * n);
422
423 for (size_t i = 0; i < n * n; i++)
424 gram [i] = gramd [i];
425
426 eigenvalues.resize (d_in);
427 // fill in only the n first ones
428 for (size_t i = 0; i < n; i++)
429 eigenvalues [i] = eigenvaluesd [i];
430
431 { // compute PCAMat = x' * v
432 FINTEGER di = d_in, ni = n;
433 float one = 1.0;
434
435 sgemm_ ("Non", "Non Trans",
436 &di, &ni, &ni,
437 &one, xc.data(), &di, gram.data(), &ni,
438 &one, PCAMat.data(), &di);
439 }
440
441 if(verbose && d_in <= 10) {
442 float *ci = PCAMat.data();
443 printf("PCAMat=\n");
444 for(int i = 0; i < n; i++) {
445 for(int j = 0; j < d_in; j++)
446 printf("%10g ", *ci++);
447 printf("\n");
448 }
449 }
450 fvec_renorm_L2 (d_in, n, PCAMat.data());
451
452 }
453
454 prepare_Ab();
455 is_trained = true;
456 }
457
458 void PCAMatrix::copy_from (const PCAMatrix & other)
459 {
460 FAISS_THROW_IF_NOT (other.is_trained);
461 mean = other.mean;
462 eigenvalues = other.eigenvalues;
463 PCAMat = other.PCAMat;
464 prepare_Ab ();
465 is_trained = true;
466 }
467
468 void PCAMatrix::prepare_Ab ()
469 {
470 FAISS_THROW_IF_NOT_FMT (
471 d_out * d_in <= PCAMat.size(),
472 "PCA matrix cannot output %d dimensions from %d ",
473 d_out, d_in);
474
475 if (!random_rotation) {
476 A = PCAMat;
477 A.resize(d_out * d_in); // strip off useless dimensions
478
479 // first scale the components
480 if (eigen_power != 0) {
481 float *ai = A.data();
482 for (int i = 0; i < d_out; i++) {
483 float factor = pow(eigenvalues[i], eigen_power);
484 for(int j = 0; j < d_in; j++)
485 *ai++ *= factor;
486 }
487 }
488
489 if (balanced_bins != 0) {
490 FAISS_THROW_IF_NOT (d_out % balanced_bins == 0);
491 int dsub = d_out / balanced_bins;
492 std::vector <float> Ain;
493 std::swap(A, Ain);
494 A.resize(d_out * d_in);
495
496 std::vector <float> accu(balanced_bins);
497 std::vector <int> counter(balanced_bins);
498
499 // greedy assignment
500 for (int i = 0; i < d_out; i++) {
501 // find best bin
502 int best_j = -1;
503 float min_w = 1e30;
504 for (int j = 0; j < balanced_bins; j++) {
505 if (counter[j] < dsub && accu[j] < min_w) {
506 min_w = accu[j];
507 best_j = j;
508 }
509 }
510 int row_dst = best_j * dsub + counter[best_j];
511 accu[best_j] += eigenvalues[i];
512 counter[best_j] ++;
513 memcpy (&A[row_dst * d_in], &Ain[i * d_in],
514 d_in * sizeof (A[0]));
515 }
516
517 if (verbose) {
518 printf(" bin accu=[");
519 for (int i = 0; i < balanced_bins; i++)
520 printf("%g ", accu[i]);
521 printf("]\n");
522 }
523 }
524
525
526 } else {
527 FAISS_THROW_IF_NOT_MSG (balanced_bins == 0,
528 "both balancing bins and applying a random rotation "
529 "does not make sense");
530 RandomRotationMatrix rr(d_out, d_out);
531
532 rr.init(5);
533
534 // apply scaling on the rotation matrix (right multiplication)
535 if (eigen_power != 0) {
536 for (int i = 0; i < d_out; i++) {
537 float factor = pow(eigenvalues[i], eigen_power);
538 for(int j = 0; j < d_out; j++)
539 rr.A[j * d_out + i] *= factor;
540 }
541 }
542
543 A.resize(d_in * d_out);
544 {
545 FINTEGER dii = d_in, doo = d_out;
546 float one = 1.0, zero = 0.0;
547
548 sgemm_ ("Not", "Not", &dii, &doo, &doo,
549 &one, PCAMat.data(), &dii, rr.A.data(), &doo, &zero,
550 A.data(), &dii);
551
552 }
553
554 }
555
556 b.clear(); b.resize(d_out);
557
558 for (int i = 0; i < d_out; i++) {
559 float accu = 0;
560 for (int j = 0; j < d_in; j++)
561 accu -= mean[j] * A[j + i * d_in];
562 b[i] = accu;
563 }
564
565 is_orthonormal = eigen_power == 0;
566
567 }
568
569 /*********************************************
570 * OPQMatrix
571 *********************************************/
572
573
574 OPQMatrix::OPQMatrix (int d, int M, int d2):
575 LinearTransform (d, d2 == -1 ? d : d2, false), M(M),
576 niter (50),
577 niter_pq (4), niter_pq_0 (40),
578 verbose(false),
579 pq(nullptr)
580 {
581 is_trained = false;
582 // OPQ is quite expensive to train, so set this right.
583 max_train_points = 256 * 256;
584 pq = nullptr;
585 }
586
587
588
589 void OPQMatrix::train (Index::idx_t n, const float *x)
590 {
591
592 const float * x_in = x;
593
594 x = fvecs_maybe_subsample (d_in, (size_t*)&n,
595 max_train_points, x, verbose);
596
597 ScopeDeleter<float> del_x (x != x_in ? x : nullptr);
598
599 // To support d_out > d_in, we pad input vectors with 0s to d_out
600 size_t d = d_out <= d_in ? d_in : d_out;
601 size_t d2 = d_out;
602
603 #if 0
604 // what this test shows: the only way of getting bit-exact
605 // reproducible results with sgeqrf and sgesvd seems to be forcing
606 // single-threading.
607 { // test repro
608 std::vector<float> r (d * d);
609 float * rotation = r.data();
610 float_randn (rotation, d * d, 1234);
611 printf("CS0: %016lx\n",
612 ivec_checksum (128*128, (int*)rotation));
613 matrix_qr (d, d, rotation);
614 printf("CS1: %016lx\n",
615 ivec_checksum (128*128, (int*)rotation));
616 return;
617 }
618 #endif
619
620 if (verbose) {
621 printf ("OPQMatrix::train: training an OPQ rotation matrix "
622 "for M=%d from %ld vectors in %dD -> %dD\n",
623 M, n, d_in, d_out);
624 }
625
626 std::vector<float> xtrain (n * d);
627 // center x
628 {
629 std::vector<float> sum (d);
630 const float *xi = x;
631 for (size_t i = 0; i < n; i++) {
632 for (int j = 0; j < d_in; j++)
633 sum [j] += *xi++;
634 }
635 for (int i = 0; i < d; i++) sum[i] /= n;
636 float *yi = xtrain.data();
637 xi = x;
638 for (size_t i = 0; i < n; i++) {
639 for (int j = 0; j < d_in; j++)
640 *yi++ = *xi++ - sum[j];
641 yi += d - d_in;
642 }
643 }
644 float *rotation;
645
646 if (A.size () == 0) {
647 A.resize (d * d);
648 rotation = A.data();
649 if (verbose)
650 printf(" OPQMatrix::train: making random %ld*%ld rotation\n",
651 d, d);
652 float_randn (rotation, d * d, 1234);
653 matrix_qr (d, d, rotation);
654 // we use only the d * d2 upper part of the matrix
655 A.resize (d * d2);
656 } else {
657 FAISS_THROW_IF_NOT (A.size() == d * d2);
658 rotation = A.data();
659 }
660
661 std::vector<float>
662 xproj (d2 * n), pq_recons (d2 * n), xxr (d * n),
663 tmp(d * d * 4);
664
665 std::vector<uint8_t> codes (M * n);
666
667 ProductQuantizer pq_default (d2, M, 8);
668 ProductQuantizer &pq_regular =
669 pq ? *pq : pq_default;
670 double t0 = getmillisecs();
671 for (int iter = 0; iter < niter; iter++) {
672
673 { // torch.mm(xtrain, rotation:t())
674 FINTEGER di = d, d2i = d2, ni = n;
675 float zero = 0, one = 1;
676 sgemm_ ("Transposed", "Not transposed",
677 &d2i, &ni, &di,
678 &one, rotation, &di,
679 xtrain.data(), &di,
680 &zero, xproj.data(), &d2i);
681 }
682
683 pq_regular.cp.max_points_per_centroid = 1000;
684 pq_regular.cp.niter = iter == 0 ? niter_pq_0 : niter_pq;
685 pq_regular.cp.verbose = verbose;
686 pq_regular.train (n, xproj.data());
687
688 pq_regular.compute_codes (xproj.data(), codes.data(), n);
689 pq_regular.decode (codes.data(), pq_recons.data(), n);
690
691 float pq_err = fvec_L2sqr (pq_recons.data(), xproj.data(), n * d2) / n;
692
693 if (verbose)
694 printf (" Iteration %d (%d PQ iterations):"
695 "%.3f s, obj=%g\n", iter, pq_regular.cp.niter,
696 (getmillisecs () - t0) / 1000.0, pq_err);
697
698 {
699 float *u = tmp.data(), *vt = &tmp [d * d];
700 float *sing_val = &tmp [2 * d * d];
701 FINTEGER di = d, d2i = d2, ni = n;
702 float one = 1, zero = 0;
703
704 // torch.mm(xtrain:t(), pq_recons)
705 sgemm_ ("Not", "Transposed",
706 &d2i, &di, &ni,
707 &one, pq_recons.data(), &d2i,
708 xtrain.data(), &di,
709 &zero, xxr.data(), &d2i);
710
711
712 FINTEGER lwork = -1, info = -1;
713 float worksz;
714 // workspace query
715 sgesvd_ ("All", "All",
716 &d2i, &di, xxr.data(), &d2i,
717 sing_val,
718 vt, &d2i, u, &di,
719 &worksz, &lwork, &info);
720
721 lwork = int(worksz);
722 std::vector<float> work (lwork);
723 // u and vt swapped
724 sgesvd_ ("All", "All",
725 &d2i, &di, xxr.data(), &d2i,
726 sing_val,
727 vt, &d2i, u, &di,
728 work.data(), &lwork, &info);
729
730 sgemm_ ("Transposed", "Transposed",
731 &di, &d2i, &d2i,
732 &one, u, &di, vt, &d2i,
733 &zero, rotation, &di);
734
735 }
736 pq_regular.train_type = ProductQuantizer::Train_hot_start;
737 }
738
739 // revert A matrix
740 if (d > d_in) {
741 for (long i = 0; i < d_out; i++)
742 memmove (&A[i * d_in], &A[i * d], sizeof(A[0]) * d_in);
743 A.resize (d_in * d_out);
744 }
745
746 is_trained = true;
747 is_orthonormal = true;
748 }
749
750
751 /*********************************************
752 * NormalizationTransform
753 *********************************************/
754
755 NormalizationTransform::NormalizationTransform (int d, float norm):
756 VectorTransform (d, d), norm (norm)
757 {
758 }
759
760 NormalizationTransform::NormalizationTransform ():
761 VectorTransform (-1, -1), norm (-1)
762 {
763 }
764
765 void NormalizationTransform::apply_noalloc
766 (idx_t n, const float* x, float* xt) const
767 {
768 if (norm == 2.0) {
769 memcpy (xt, x, sizeof (x[0]) * n * d_in);
770 fvec_renorm_L2 (d_in, n, xt);
771 } else {
772 FAISS_THROW_MSG ("not implemented");
773 }
774 }
775
776 void NormalizationTransform::reverse_transform (idx_t n, const float* xt,
777 float* x) const
778 {
779 memcpy (x, xt, sizeof (xt[0]) * n * d_in);
780 }
781
782 /*********************************************
783 * IndexPreTransform
784 *********************************************/
785
786 IndexPreTransform::IndexPreTransform ():
787 index(nullptr), own_fields (false)
788 {
789 }
790
791
792 IndexPreTransform::IndexPreTransform (
793 Index * index):
794 Index (index->d, index->metric_type),
795 index (index), own_fields (false)
796 {
797 is_trained = index->is_trained;
798 }
799
800
801 IndexPreTransform::IndexPreTransform (
802 VectorTransform * ltrans,
803 Index * index):
804 Index (index->d, index->metric_type),
805 index (index), own_fields (false)
806 {
807 is_trained = index->is_trained;
808 prepend_transform (ltrans);
809 }
810
811 void IndexPreTransform::prepend_transform (VectorTransform *ltrans)
812 {
813 FAISS_THROW_IF_NOT (ltrans->d_out == d);
814 is_trained = is_trained && ltrans->is_trained;
815 chain.insert (chain.begin(), ltrans);
816 d = ltrans->d_in;
817 }
818
819
820 IndexPreTransform::~IndexPreTransform ()
821 {
822 if (own_fields) {
823 for (int i = 0; i < chain.size(); i++)
824 delete chain[i];
825 delete index;
826 }
827 }
828
829
830
831
832 void IndexPreTransform::train (idx_t n, const float *x)
833 {
834 int last_untrained = 0;
835 if (!index->is_trained) {
836 last_untrained = chain.size();
837 } else {
838 for (int i = chain.size() - 1; i >= 0; i--) {
839 if (!chain[i]->is_trained) {
840 last_untrained = i;
841 break;
842 }
843 }
844 }
845 const float *prev_x = x;
846 ScopeDeleter<float> del;
847
848 if (verbose) {
849 printf("IndexPreTransform::train: training chain 0 to %d\n",
850 last_untrained);
851 }
852
853 for (int i = 0; i <= last_untrained; i++) {
854
855 if (i < chain.size()) {
856 VectorTransform *ltrans = chain [i];
857 if (!ltrans->is_trained) {
858 if (verbose) {
859 printf(" Training chain component %d/%zd\n",
860 i, chain.size());
861 if (OPQMatrix *opqm = dynamic_cast<OPQMatrix*>(ltrans)) {
862 opqm->verbose = true;
863 }
864 }
865 ltrans->train (n, prev_x);
866 }
867 } else {
868 if (verbose) {
869 printf(" Training sub-index\n");
870 }
871 index->train (n, prev_x);
872 }
873 if (i == last_untrained) break;
874 if (verbose) {
875 printf(" Applying transform %d/%zd\n",
876 i, chain.size());
877 }
878
879 float * xt = chain[i]->apply (n, prev_x);
880
881 if (prev_x != x) delete [] prev_x;
882 prev_x = xt;
883 del.set(xt);
884 }
885
886 is_trained = true;
887 }
888
889
890 const float *IndexPreTransform::apply_chain (idx_t n, const float *x) const
891 {
892 const float *prev_x = x;
893 ScopeDeleter<float> del;
894
895 for (int i = 0; i < chain.size(); i++) {
896 float * xt = chain[i]->apply (n, prev_x);
897 ScopeDeleter<float> del2 (xt);
898 del2.swap (del);
899 prev_x = xt;
900 }
901 del.release ();
902 return prev_x;
903 }
904
905 void IndexPreTransform::reverse_chain (idx_t n, const float* xt, float* x) const
906 {
907 const float* next_x = xt;
908 ScopeDeleter<float> del;
909
910 for (int i = chain.size() - 1; i >= 0; i--) {
911 float* prev_x = (i == 0) ? x : new float [n * chain[i]->d_in];
912 ScopeDeleter<float> del2 ((prev_x == x) ? nullptr : prev_x);
913 chain [i]->reverse_transform (n, next_x, prev_x);
914 del2.swap (del);
915 next_x = prev_x;
916 }
917 }
918
919 void IndexPreTransform::add (idx_t n, const float *x)
920 {
921 FAISS_THROW_IF_NOT (is_trained);
922 const float *xt = apply_chain (n, x);
923 ScopeDeleter<float> del(xt == x ? nullptr : xt);
924 index->add (n, xt);
925 ntotal = index->ntotal;
926 }
927
928 void IndexPreTransform::add_with_ids (idx_t n, const float * x,
929 const long *xids)
930 {
931 FAISS_THROW_IF_NOT (is_trained);
932 const float *xt = apply_chain (n, x);
933 ScopeDeleter<float> del(xt == x ? nullptr : xt);
934 index->add_with_ids (n, xt, xids);
935 ntotal = index->ntotal;
936 }
937
938
939
940
941 void IndexPreTransform::search (idx_t n, const float *x, idx_t k,
942 float *distances, idx_t *labels) const
943 {
944 FAISS_THROW_IF_NOT (is_trained);
945 const float *xt = apply_chain (n, x);
946 ScopeDeleter<float> del(xt == x ? nullptr : xt);
947 index->search (n, xt, k, distances, labels);
948 }
949
950
951 void IndexPreTransform::reset () {
952 index->reset();
953 ntotal = 0;
954 }
955
956 long IndexPreTransform::remove_ids (const IDSelector & sel) {
957 long nremove = index->remove_ids (sel);
958 ntotal = index->ntotal;
959 return nremove;
960 }
961
962
963 void IndexPreTransform::reconstruct (idx_t key, float * recons) const
964 {
965 float *x = chain.empty() ? recons : new float [index->d];
966 ScopeDeleter<float> del (recons == x ? nullptr : x);
967 // Initial reconstruction
968 index->reconstruct (key, x);
969
970 // Revert transformations from last to first
971 reverse_chain (1, x, recons);
972 }
973
974
975 void IndexPreTransform::reconstruct_n (idx_t i0, idx_t ni, float *recons) const
976 {
977 float *x = chain.empty() ? recons : new float [ni * index->d];
978 ScopeDeleter<float> del (recons == x ? nullptr : x);
979 // Initial reconstruction
980 index->reconstruct_n (i0, ni, x);
981
982 // Revert transformations from last to first
983 reverse_chain (ni, x, recons);
984 }
985
986
987 void IndexPreTransform::search_and_reconstruct (
988 idx_t n, const float *x, idx_t k,
989 float *distances, idx_t *labels, float* recons) const
990 {
991 FAISS_THROW_IF_NOT (is_trained);
992
993 const float* xt = apply_chain (n, x);
994 ScopeDeleter<float> del ((xt == x) ? nullptr : xt);
995
996 float* recons_temp = chain.empty() ? recons : new float [n * k * index->d];
997 ScopeDeleter<float> del2 ((recons_temp == recons) ? nullptr : recons_temp);
998 index->search_and_reconstruct (n, xt, k, distances, labels, recons_temp);
999
1000 // Revert transformations from last to first
1001 reverse_chain (n * k, recons_temp, recons);
1002 }
1003
1004
1005 /*********************************************
1006 * RemapDimensionsTransform
1007 *********************************************/
1008
1009
1010 RemapDimensionsTransform::RemapDimensionsTransform (
1011 int d_in, int d_out, const int *map_in):
1012 VectorTransform (d_in, d_out)
1013 {
1014 map.resize (d_out);
1015 for (int i = 0; i < d_out; i++) {
1016 map[i] = map_in[i];
1017 FAISS_THROW_IF_NOT (map[i] == -1 || (map[i] >= 0 && map[i] < d_in));
1018 }
1019 }
1020
1021 RemapDimensionsTransform::RemapDimensionsTransform (
1022 int d_in, int d_out, bool uniform): VectorTransform (d_in, d_out)
1023 {
1024 map.resize (d_out, -1);
1025
1026 if (uniform) {
1027 if (d_in < d_out) {
1028 for (int i = 0; i < d_in; i++) {
1029 map [i * d_out / d_in] = i;
1030 }
1031 } else {
1032 for (int i = 0; i < d_out; i++) {
1033 map [i] = i * d_in / d_out;
1034 }
1035 }
1036 } else {
1037 for (int i = 0; i < d_in && i < d_out; i++)
1038 map [i] = i;
1039 }
1040 }
1041
1042
1043 void RemapDimensionsTransform::apply_noalloc (idx_t n, const float * x,
1044 float *xt) const
1045 {
1046 for (idx_t i = 0; i < n; i++) {
1047 for (int j = 0; j < d_out; j++) {
1048 xt[j] = map[j] < 0 ? 0 : x[map[j]];
1049 }
1050 x += d_in;
1051 xt += d_out;
1052 }
1053 }
1054
1055 void RemapDimensionsTransform::reverse_transform (idx_t n, const float * xt,
1056 float *x) const
1057 {
1058 memset (x, 0, sizeof (*x) * n * d_in);
1059 for (idx_t i = 0; i < n; i++) {
1060 for (int j = 0; j < d_out; j++) {
1061 if (map[j] >= 0) x[map[j]] = xt[j];
1062 }
1063 x += d_in;
1064 xt += d_out;
1065 }
1066 }
faiss::LinearTransform::transform_transpose
void transform_transpose(idx_t n, const float *y, float *x) const
Definition: VectorTransform.cpp:135
faiss::IndexPreTransform::index
Index * index
! chain of tranforms
Definition: VectorTransform.h:254
faiss::RandomRotationMatrix
Randomly rotate a set of vectors.
Definition: VectorTransform.h:117
faiss::ClusteringParameters::niter
int niter
clustering iterations
Definition: Clustering.h:25
faiss::OPQMatrix::niter
int niter
Number of outer training iterations.
Definition: VectorTransform.h:190
faiss::ProductQuantizer::decode
void decode(const uint8_t *code, float *x) const
decode a vector from a given code (or n vectors if third argument)
Definition: ProductQuantizer.cpp:339
faiss::fvec_L2sqr
float fvec_L2sqr(const float *x, const float *y, size_t d)
Squared L2 distance between two vectors.
Definition: utils.cpp:574
faiss::RandomRotationMatrix::init
void init(int seed)
must be called before the transform is used
Definition: VectorTransform.cpp:215
faiss::IndexPreTransform::reset
void reset() override
removes all elements from the database.
Definition: VectorTransform.cpp:951
faiss::Index::reset
virtual void reset()=0
removes all elements from the database.
faiss::OPQMatrix::niter_pq
int niter_pq
Number of training iterations for the PQ.
Definition: VectorTransform.h:191
faiss::LinearTransform::A
std::vector< float > A
Transformation matrix, size d_out * d_in.
Definition: VectorTransform.h:85
faiss::Index::search_and_reconstruct
virtual void search_and_reconstruct(idx_t n, const float *x, idx_t k, float *distances, idx_t *labels, float *recons) const
Definition: Index.cpp:66
faiss::fvecs_maybe_subsample
const float * fvecs_maybe_subsample(size_t d, size_t *n, size_t nmax, const float *x, bool verbose, long seed)
Definition: utils.cpp:1975
faiss::PCAMatrix::eigen_power
float eigen_power
Definition: VectorTransform.h:140
faiss::IDSelector
Definition: AuxIndexStructures.h:52
faiss::LinearTransform::LinearTransform
LinearTransform(int d_in=0, int d_out=0, bool have_bias=false)
both d_in &gt; d_out and d_out &lt; d_in are supported
Definition: VectorTransform.cpp:100
faiss::PCAMatrix
Definition: VectorTransform.h:132
faiss::LinearTransform::set_is_orthonormal
void set_is_orthonormal()
compute A^T * A to set the is_orthonormal flag
Definition: VectorTransform.cpp:160
faiss::OPQMatrix::pq
ProductQuantizer * pq
Definition: VectorTransform.h:200
faiss::Index::train
virtual void train(idx_t n, const float *x)
Definition: Index.cpp:23
faiss::Index::add_with_ids
virtual void add_with_ids(idx_t n, const float *x, const long *xids)
Definition: Index.cpp:41
faiss::OPQMatrix
Definition: VectorTransform.h:187
faiss::ScopeDeleter
Definition: FaissException.h:39
faiss::PCAMatrix::train
void train(Index::idx_t n, const float *x) override
Definition: VectorTransform.cpp:316
faiss::PCAMatrix::mean
std::vector< float > mean
Mean, size d_in.
Definition: VectorTransform.h:152
faiss::RemapDimensionsTransform::map
std::vector< int > map
Definition: VectorTransform.h:216
faiss::IndexPreTransform::apply_chain
const float * apply_chain(idx_t n, const float *x) const
Definition: VectorTransform.cpp:890
faiss::PCAMatrix::PCAMat
std::vector< float > PCAMat
PCA matrix, size d_in * d_in.
Definition: VectorTransform.h:158
faiss::IndexPreTransform::train
void train(idx_t n, const float *x) override
Definition: VectorTransform.cpp:832
faiss::ProductQuantizer::compute_codes
void compute_codes(const float *x, uint8_t *codes, size_t n) const
same as compute_code for several vectors
Definition: ProductQuantizer.cpp:386
faiss::Index::d
int d
vector dimension
Definition: Index.h:64
faiss::LinearTransform::b
std::vector< float > b
bias vector, size d_out
Definition: VectorTransform.h:88
faiss::LinearTransform::reverse_transform
void reverse_transform(idx_t n, const float *xt, float *x) const override
works only if is_orthonormal
Definition: VectorTransform.cpp:199
faiss::Index::reconstruct_n
virtual void reconstruct_n(idx_t i0, idx_t ni, float *recons) const
Definition: Index.cpp:59
faiss::Index::add
virtual void add(idx_t n, const float *x)=0
faiss::NormalizationTransform::reverse_transform
void reverse_transform(idx_t n, const float *xt, float *x) const override
Identity transform since norm is not revertible.
Definition: VectorTransform.cpp:776
faiss::LinearTransform
Definition: VectorTransform.h:77
faiss::OPQMatrix::train
void train(Index::idx_t n, const float *x) override
Definition: VectorTransform.cpp:589
faiss::PCAMatrix::balanced_bins
int balanced_bins
try to distribute output eigenvectors in this many bins
Definition: VectorTransform.h:149
faiss::VectorTransform::is_trained
bool is_trained
Definition: VectorTransform.h:42
faiss::Index::idx_t
long idx_t
all indices are this type
Definition: Index.h:62
faiss::IndexPreTransform::reconstruct_n
void reconstruct_n(idx_t i0, idx_t ni, float *recons) const override
Definition: VectorTransform.cpp:975
faiss::Index::ntotal
idx_t ntotal
total nb of indexed vectors
Definition: Index.h:65
faiss::NormalizationTransform::apply_noalloc
void apply_noalloc(idx_t n, const float *x, float *xt) const override
same as apply, but result is pre-allocated
Definition: VectorTransform.cpp:766
faiss::ProductQuantizer::Train_hot_start
the centroids are already initialized
Definition: ProductQuantizer.h:42
faiss::getmillisecs
double getmillisecs()
ms elapsed since some arbitrary epoch
Definition: utils.cpp:74
faiss::Index::remove_ids
virtual long remove_ids(const IDSelector &sel)
Definition: Index.cpp:48
faiss::Index::search
virtual void search(idx_t n, const float *x, idx_t k, float *distances, idx_t *labels) const =0
faiss::VectorTransform
Definition: VectorTransform.h:28
faiss::matrix_qr
void matrix_qr(int m, int n, float *a)
Definition: utils.cpp:1322
faiss::IndexPreTransform::own_fields
bool own_fields
! the sub-index
Definition: VectorTransform.h:256
faiss::OPQMatrix::niter_pq_0
int niter_pq_0
same, for the first outer iteration
Definition: VectorTransform.h:192
faiss::ProductQuantizer::cp
ClusteringParameters cp
parameters used during clustering
Definition: ProductQuantizer.h:49
faiss::ivec_checksum
size_t ivec_checksum(size_t n, const int *a)
compute a checksum on a table.
Definition: utils.cpp:1672
faiss::VectorTransform::train
virtual void train(idx_t n, const float *x)
Definition: VectorTransform.cpp:80
faiss::VectorTransform::reverse_transform
virtual void reverse_transform(idx_t n, const float *xt, float *x) const
Definition: VectorTransform.cpp:85
faiss::IndexPreTransform::search_and_reconstruct
void search_and_reconstruct(idx_t n, const float *x, idx_t k, float *distances, idx_t *labels, float *recons) const override
Definition: VectorTransform.cpp:987
faiss::RemapDimensionsTransform::reverse_transform
void reverse_transform(idx_t n, const float *xt, float *x) const override
reverse transform correct only when the mapping is a permuation
Definition: VectorTransform.cpp:1055
faiss::OPQMatrix::max_train_points
size_t max_train_points
if there are too many training points, resample
Definition: VectorTransform.h:195
faiss::PCAMatrix::copy_from
void copy_from(const PCAMatrix &other)
copy pre-trained PCA matrix
Definition: VectorTransform.cpp:458
faiss::VectorTransform::d_out
int d_out
! input dimension
Definition: VectorTransform.h:33
faiss::OPQMatrix::OPQMatrix
OPQMatrix(int d=0, int M=1, int d2=-1)
if d2 != -1, output vectors of this dimension
Definition: VectorTransform.cpp:574
faiss::Index
Definition: Index.h:60
faiss::PCAMatrix::prepare_Ab
void prepare_Ab()
called after mean, PCAMat and eigenvalues are computed
Definition: VectorTransform.cpp:468
faiss::IndexPreTransform::add
void add(idx_t n, const float *x) override
Definition: VectorTransform.cpp:919
faiss::RemapDimensionsTransform::apply_noalloc
void apply_noalloc(idx_t n, const float *x, float *xt) const override
same as apply, but result is pre-allocated
Definition: VectorTransform.cpp:1043
faiss::ProductQuantizer
Definition: ProductQuantizer.h:25
faiss::Index::is_trained
bool is_trained
set if the Index does not require training, or if training is done already
Definition: Index.h:69
faiss::IndexPreTransform::reverse_chain
void reverse_chain(idx_t n, const float *xt, float *x) const
Definition: VectorTransform.cpp:905
faiss::LinearTransform::is_orthonormal
bool is_orthonormal
! whether to use the bias term
Definition: VectorTransform.h:82
faiss::PCAMatrix::eigenvalues
std::vector< float > eigenvalues
eigenvalues of covariance matrix (= squared singular values)
Definition: VectorTransform.h:155
faiss::IndexPreTransform::search
void search(idx_t n, const float *x, idx_t k, float *distances, idx_t *labels) const override
Definition: VectorTransform.cpp:941
faiss::Index::reconstruct
virtual void reconstruct(idx_t key, float *recons) const
Definition: Index.cpp:54
faiss::IndexPreTransform::add_with_ids
void add_with_ids(idx_t n, const float *x, const long *xids) override
Definition: VectorTransform.cpp:928
faiss::PCAMatrix::random_rotation
bool random_rotation
random rotation after PCA
Definition: VectorTransform.h:143
faiss::PCAMatrix::max_points_per_d
size_t max_points_per_d
ratio between # training vectors and dimension
Definition: VectorTransform.h:146
faiss::ClusteringParameters::max_points_per_centroid
int max_points_per_centroid
to limit size of dataset
Definition: Clustering.h:34
faiss::VectorTransform::apply
float * apply(idx_t n, const float *x) const
Definition: VectorTransform.cpp:72
faiss::IndexPreTransform::remove_ids
long remove_ids(const IDSelector &sel) override
Definition: VectorTransform.cpp:956
faiss::VectorTransform::apply_noalloc
virtual void apply_noalloc(idx_t n, const float *x, float *xt) const =0
same as apply, but result is pre-allocated
faiss::IndexPreTransform::reconstruct
void reconstruct(idx_t key, float *recons) const override
Definition: VectorTransform.cpp:963
faiss::OPQMatrix::M
int M
nb of subquantizers
Definition: VectorTransform.h:189
faiss::LinearTransform::apply_noalloc
void apply_noalloc(idx_t n, const float *x, float *xt) const override
same as apply, but result is pre-allocated
Definition: VectorTransform.cpp:106