Faiss

- Main Page
- Namespaces
- Classes
- Files

- File List

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