001/**
002 * Copyright (c) 2011, The University of Southampton and the individual contributors.
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without modification,
006 * are permitted provided that the following conditions are met:
007 *
008 *   *  Redistributions of source code must retain the above copyright notice,
009 *      this list of conditions and the following disclaimer.
010 *
011 *   *  Redistributions in binary form must reproduce the above copyright notice,
012 *      this list of conditions and the following disclaimer in the documentation
013 *      and/or other materials provided with the distribution.
014 *
015 *   *  Neither the name of the University of Southampton nor the names of its
016 *      contributors may be used to endorse or promote products derived from this
017 *      software without specific prior written permission.
018 *
019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
029 */
030package org.openimaj.math.matrix.algorithm.pca;
031
032import java.util.Arrays;
033import java.util.List;
034
035import Jama.Matrix;
036
037/**
038 * Abstract base class for PCA implementations.
039 * 
040 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
041 * 
042 */
043public abstract class PrincipalComponentAnalysis {
044        /**
045         * Interface for classes capable of selecting a subset of the PCA components
046         * 
047         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
048         */
049        public interface ComponentSelector {
050                /**
051                 * Select a subset of principal components, discarding the ones not
052                 * selected
053                 * 
054                 * @param pca
055                 */
056                public void select(PrincipalComponentAnalysis pca);
057        }
058
059        /**
060         * {@link ComponentSelector} that selects the n-best components.
061         * 
062         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
063         * 
064         */
065        public static class NumberComponentSelector implements ComponentSelector {
066                int n;
067
068                /**
069                 * Construct with the number of components required
070                 * 
071                 * @param n
072                 *            the number of components
073                 */
074                public NumberComponentSelector(int n) {
075                        this.n = n;
076                }
077
078                @Override
079                public void select(PrincipalComponentAnalysis pca) {
080                        pca.selectSubset(n);
081                }
082        }
083
084        /**
085         * {@link ComponentSelector} that selects a subset of the principal
086         * components such that all remaining components have a cumulative energy
087         * less than the given value.
088         * 
089         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
090         */
091        public static class EnergyThresholdComponentSelector implements ComponentSelector {
092                double threshold;
093
094                /**
095                 * Construct with the given threshold
096                 * 
097                 * @param threshold
098                 *            the threshold
099                 */
100                public EnergyThresholdComponentSelector(double threshold) {
101                        this.threshold = threshold;
102                }
103
104                @Override
105                public void select(PrincipalComponentAnalysis pca) {
106                        pca.selectSubsetEnergyThreshold(threshold);
107                }
108        }
109
110        /**
111         * {@link ComponentSelector} that selects a subset of the principal
112         * components such that all remaining components have a certain percentage
113         * cumulative energy of the total. The percentage is calculated relative to
114         * the total energy of the eigenvalues.
115         * 
116         * Bear in mind that if not all the eigenvalues were calculated, or if some
117         * have previously been removed through {@link #selectSubset(int)},
118         * {@link #selectSubsetEnergyThreshold(double)} or
119         * {@link #selectSubsetPercentageEnergy(double)}, then the percentage
120         * calculation only factors in the remaining eigenvalues that are available
121         * to it.
122         * 
123         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
124         */
125        public static class PercentageEnergyComponentSelector implements ComponentSelector {
126                double percentage;
127
128                /**
129                 * Construct with the given percentage
130                 * 
131                 * @param percentage
132                 *            percentage of the total cumulative energy to retain
133                 *            [0..1].
134                 */
135                public PercentageEnergyComponentSelector(double percentage) {
136                        this.percentage = percentage;
137                }
138
139                @Override
140                public void select(PrincipalComponentAnalysis pca) {
141                        pca.selectSubsetPercentageEnergy(percentage);
142                }
143        }
144
145        protected Matrix basis;
146        protected double[] mean;
147        protected double[] eigenvalues;
148
149        /**
150         * Learn the principal components of the given data list. Each item
151         * corresponds to an observation with the number of dimensions equal to the
152         * length of the array.
153         * 
154         * @param data
155         *            the data
156         */
157        public void learnBasis(List<double[]> data) {
158                learnBasis(data.toArray(new double[data.size()][]));
159        }
160
161        /**
162         * Learn the principal components of the given data array. Each row
163         * corresponds to an observation with the number of dimensions equal to the
164         * number of columns.
165         * 
166         * @param data
167         *            the data
168         */
169        public void learnBasis(double[][] data) {
170                learnBasis(new Matrix(data));
171        }
172
173        /**
174         * Learn the principal components of the given data matrix. Each row
175         * corresponds to an observation with the number of dimensions equal to the
176         * number of columns.
177         * 
178         * @param data
179         *            the data
180         */
181        public void learnBasis(Matrix data) {
182                final Matrix norm = this.buildNormalisedDataMatrix(data);
183                learnBasisNorm(norm);
184        }
185
186        /**
187         * Learn the PCA basis from the centered data provided. Each row corresponds
188         * to an observation with the number of dimensions equal to the number of
189         * columns.
190         * 
191         * @param norm
192         */
193        protected abstract void learnBasisNorm(Matrix norm);
194
195        /**
196         * Zero-centre the data matrix and return a copy
197         * 
198         * @param data
199         *            the data matrix
200         * @return the normalised data
201         */
202        protected Matrix buildNormalisedDataMatrix(Matrix m) {
203                final double[][] data = m.getArray();
204
205                mean = new double[data[0].length];
206
207                for (int j = 0; j < data.length; j++)
208                        for (int i = 0; i < data[0].length; i++)
209                                mean[i] += data[j][i];
210
211                for (int i = 0; i < data[0].length; i++)
212                        mean[i] /= data.length;
213
214                final Matrix mat = new Matrix(data.length, data[0].length);
215                final double[][] matdat = mat.getArray();
216
217                for (int j = 0; j < data.length; j++)
218                        for (int i = 0; i < data[0].length; i++)
219                                matdat[j][i] = (data[j][i] - mean[i]);
220
221                return mat;
222        }
223
224        /**
225         * Get the principal components. The principal components are the column
226         * vectors of the returned matrix.
227         * 
228         * @return the principal components
229         */
230        public Matrix getBasis() {
231                return basis;
232        }
233
234        /**
235         * Get a specific principle component vector as a double array. The returned
236         * array contains a copy of the data.
237         * 
238         * @param index
239         *            the index of the principle component
240         * 
241         * @return the principle component
242         */
243        public double[] getPrincipalComponent(int index) {
244                final double[] pc = new double[basis.getRowDimension()];
245                final double[][] data = basis.getArray();
246
247                for (int r = 0; r < pc.length; r++)
248                        pc[r] = data[r][index];
249
250                return pc;
251        }
252
253        /**
254         * Get the principal components. The principal components are the column
255         * vectors of the returned matrix.
256         * 
257         * Syntactic sugar for {@link #getBasis()}
258         * 
259         * @return the principal components
260         */
261        public Matrix getEigenVectors() {
262                return basis;
263        }
264
265        /**
266         * @return the eigen values corresponding to the principal components
267         */
268        public double[] getEigenValues() {
269                return eigenvalues;
270        }
271
272        /**
273         * Get the eigen value corresponding to the ith principal component.
274         * 
275         * @param i
276         *            the index of the component
277         * @return the eigen value corresponding to the principal component
278         */
279        public double getEigenValue(int i) {
280                return eigenvalues[i];
281        }
282
283        /**
284         * Get the cumulative energy of the ith principal component
285         * 
286         * @param i
287         *            the index of the component
288         * @return the cumulative energy of the component
289         */
290        public double getCumulativeEnergy(int i) {
291                double energy = 0;
292
293                for (int j = 0; j <= i; j++)
294                        energy += eigenvalues[i];
295
296                return energy;
297        }
298
299        /**
300         * Get the cumulative energies of each principal component
301         * 
302         * @return the cumulative energies of the components
303         */
304        public double[] getCumulativeEnergies() {
305                final double[] energy = new double[eigenvalues.length];
306
307                energy[0] = eigenvalues[0];
308                for (int j = 1; j < energy.length; j++)
309                        energy[j] += energy[j - 1] + eigenvalues[j];
310
311                return energy;
312        }
313
314        /**
315         * @return The mean values
316         */
317        public double[] getMean() {
318                return mean;
319        }
320
321        /**
322         * Select a subset of the principal components using a
323         * {@link ComponentSelector}. Calling this method throws away any extra
324         * basis vectors and eigenvalues.
325         * 
326         * @param selector
327         *            the {@link ComponentSelector} to apply
328         */
329        public void selectSubset(ComponentSelector selector) {
330                selector.select(this);
331        }
332
333        /**
334         * Select a subset of the principal components. Calling this method throws
335         * away any extra basis vectors and eigenvalues.
336         * 
337         * @param n
338         */
339        public void selectSubset(int n) {
340                if (n >= eigenvalues.length)
341                        return;
342
343                basis = basis.getMatrix(0, basis.getRowDimension() - 1, 0, n - 1);
344                eigenvalues = Arrays.copyOf(eigenvalues, n);
345        }
346
347        /**
348         * Select a subset of the principal components such that all remaining
349         * components have a cumulative energy less than the given value.
350         * 
351         * Calling this method throws away any extra basis vectors and eigenvalues.
352         * 
353         * @param threshold
354         *            threshold on the cumulative energy.
355         */
356        public void selectSubsetEnergyThreshold(double threshold) {
357                final double[] energy = getCumulativeEnergies();
358
359                for (int i = 1; i < energy.length; i++) {
360                        if (energy[i] < threshold) {
361                                selectSubset(i - 1);
362                                return;
363                        }
364                }
365        }
366
367        /**
368         * Select a subset of the principal components such that all remaining
369         * components have a certain percentage cumulative energy of the total. The
370         * percentage is calculated relative to the total energy of the eigenvalues.
371         * Bear in mind that if not all the eigenvalues were calculated, or if some
372         * have previously been removed through {@link #selectSubset(int)},
373         * {@link #selectSubsetEnergyThreshold(double)} or
374         * {@link #selectSubsetPercentageEnergy(double)}, then the percentage
375         * calculation only factors in the remaining eigenvalues that are available
376         * to it.
377         * 
378         * Calling this method throws away any extra basis vectors and eigenvalues.
379         * 
380         * @param percentage
381         *            percentage of the total cumulative energy to retain [0..1].
382         */
383        public void selectSubsetPercentageEnergy(double percentage) {
384                final double[] energy = getCumulativeEnergies();
385                final double total = energy[energy.length - 1];
386
387                for (int i = 1; i < energy.length; i++) {
388                        if (energy[i] / total > percentage) {
389                                selectSubset(i - 1);
390                                return;
391                        }
392                }
393        }
394
395        /**
396         * Generate a new "observation" as a linear combination of the principal
397         * components (PC): mean + PC * scaling.
398         * 
399         * If the scaling vector is shorter than the number of components, it will
400         * be zero-padded. If it is longer, it will be truncated.
401         * 
402         * @param scalings
403         *            the weighting for each PC
404         * @return generated observation
405         */
406        public double[] generate(double[] scalings) {
407                final Matrix scale = new Matrix(this.eigenvalues.length, 1);
408
409                for (int i = 0; i < Math.min(eigenvalues.length, scalings.length); i++)
410                        scale.set(i, 0, scalings[i]);
411
412                final Matrix meanMatrix = new Matrix(new double[][] { mean }).transpose();
413
414                return meanMatrix.plus(basis.times(scale)).getColumnPackedCopy();
415        }
416
417        /**
418         * Project a matrix of row vectors by the basis. The vectors are normalised
419         * by subtracting the mean and then multiplied by the basis. The returned
420         * matrix has a row for each vector.
421         * 
422         * @param m
423         *            the vector to project
424         * @return projected vectors
425         */
426        public Matrix project(Matrix m) {
427                final Matrix vec = m.copy();
428
429                final int rows = vec.getRowDimension();
430                final int cols = vec.getColumnDimension();
431                final double[][] vecarr = vec.getArray();
432
433                for (int r = 0; r < rows; r++)
434                        for (int c = 0; c < cols; c++)
435                                vecarr[r][c] -= mean[c];
436
437                // T = (Vt.Dt)^T == Dt.Vt
438                return vec.times(basis);
439        }
440
441        /**
442         * Project a vector by the basis. The vector is normalised by subtracting
443         * the mean and then multiplied by the basis.
444         * 
445         * @param vector
446         *            the vector to project
447         * @return projected vector
448         */
449        public double[] project(double[] vector) {
450                final Matrix vec = new Matrix(1, vector.length);
451                final double[][] vecarr = vec.getArray();
452
453                for (int i = 0; i < vector.length; i++)
454                        vecarr[0][i] = vector[i] - mean[i];
455
456                return vec.times(basis).getColumnPackedCopy();
457        }
458
459        /**
460         * Get the standard deviations (sqrt of eigenvalues) of the principal
461         * components.
462         * 
463         * @return vector of standard deviations
464         */
465        public double[] getStandardDeviations() {
466                return getStandardDeviations(eigenvalues.length);
467        }
468
469        /**
470         * Get the standard deviations (sqrt of eigenvalues) of the n top principal
471         * components.
472         * 
473         * @param n
474         *            number of principal components
475         * @return vector of standard deviations
476         */
477        public double[] getStandardDeviations(int n) {
478                final double[] rngs = new double[n];
479
480                for (int i = 0; i < rngs.length; i++) {
481                        rngs[i] = Math.sqrt(eigenvalues[i]);
482                }
483
484                return rngs;
485        }
486
487        @Override
488        public String toString() {
489                return String.format("PrincipalComponentAnalysis[dims=%d]", this.eigenvalues == null ? 0
490                                : this.eigenvalues.length);
491        }
492}