1 /**
2 * @file methods/radical/radical.cpp
3 * @author Nishant Mehta
4 *
5 * Implementation of Radical class
6 *
7 * mlpack is free software; you may redistribute it and/or modify it under the
8 * terms of the 3-clause BSD license. You should have received a copy of the
9 * 3-clause BSD license along with mlpack. If not, see
10 * http://www.opensource.org/licenses/BSD-3-Clause for more information.
11 */
12
13 #include "radical.hpp"
14 #include <mlpack/core/util/log.hpp>
15 #include <mlpack/core/util/timers.hpp>
16
17 using namespace std;
18 using namespace arma;
19 using namespace mlpack;
20 using namespace mlpack::radical;
21
22 // Set the parameters to RADICAL.
Radical(const double noiseStdDev,const size_t replicates,const size_t angles,const size_t sweeps,const size_t m)23 Radical::Radical(const double noiseStdDev,
24 const size_t replicates,
25 const size_t angles,
26 const size_t sweeps,
27 const size_t m) :
28 noiseStdDev(noiseStdDev),
29 replicates(replicates),
30 angles(angles),
31 sweeps(sweeps),
32 m(m)
33 {
34 // Nothing to do here.
35 }
36
CopyAndPerturb(mat & xNew,const mat & x) const37 void Radical::CopyAndPerturb(mat& xNew, const mat& x) const
38 {
39 Timer::Start("radical_copy_and_perturb");
40 xNew = repmat(x, replicates, 1) + noiseStdDev * randn(replicates * x.n_rows,
41 x.n_cols);
42 Timer::Stop("radical_copy_and_perturb");
43 }
44
45
Vasicek(vec & z) const46 double Radical::Vasicek(vec& z) const
47 {
48 z = sort(z);
49
50 // Apparently slower.
51 /*
52 vec logs = log(z.subvec(m, z.n_elem - 1) - z.subvec(0, z.n_elem - 1 - m));
53 //vec val = sum(log(z.subvec(m, nPoints - 1) - z.subvec(0, nPoints - 1 - m)));
54 return (double) sum(logs);
55 */
56
57 // Apparently faster.
58 double sum = 0;
59 uword range = z.n_elem - m;
60 for (uword i = 0; i < range; ++i)
61 {
62 sum += log(max(z(i + m) - z(i), DBL_MIN));
63 }
64
65 return sum;
66 }
67
68
DoRadical2D(const mat & matX)69 double Radical::DoRadical2D(const mat& matX)
70 {
71 CopyAndPerturb(perturbed, matX);
72
73 mat::fixed<2, 2> matJacobi;
74
75 vec values(angles);
76
77 for (size_t i = 0; i < angles; ++i)
78 {
79 const double theta = (i / (double) angles) * M_PI / 2.0;
80 const double cosTheta = cos(theta);
81 const double sinTheta = sin(theta);
82
83 matJacobi(0, 0) = cosTheta;
84 matJacobi(1, 0) = -sinTheta;
85 matJacobi(0, 1) = sinTheta;
86 matJacobi(1, 1) = cosTheta;
87
88 candidate = perturbed * matJacobi;
89 vec candidateY1 = candidate.unsafe_col(0);
90 vec candidateY2 = candidate.unsafe_col(1);
91
92 values(i) = Vasicek(candidateY1) + Vasicek(candidateY2);
93 }
94
95 uword indOpt = 0;
96 values.min(indOpt); // we ignore the return value; we don't care about it
97 return (indOpt / (double) angles) * M_PI / 2.0;
98 }
99
100
DoRadical(const mat & matXT,mat & matY,mat & matW)101 void Radical::DoRadical(const mat& matXT, mat& matY, mat& matW)
102 {
103 // matX is nPoints by nDims (although less intuitive than columns being
104 // points, and although this is the transpose of the ICA literature, this
105 // choice is for computational efficiency when repeatedly generating
106 // two-dimensional coordinate projections for Radical2D).
107 Timer::Start("radical_transpose_data");
108 mat matX = trans(matXT);
109 Timer::Stop("radical_transpose_data");
110
111 // If m was not specified, initialize m as recommended in
112 // (Learned-Miller and Fisher, 2003).
113 if (m < 1)
114 m = floor(sqrt((double) matX.n_rows));
115
116 const size_t nDims = matX.n_cols;
117 const size_t nPoints = matX.n_rows;
118
119 Timer::Start("radical_whiten_data");
120 mat matXWhitened;
121 mat matWhitening;
122 WhitenFeatureMajorMatrix(matX, matY, matWhitening);
123 Timer::Stop("radical_whiten_data");
124 // matY is now the whitened form of matX.
125
126 // In the RADICAL code, they do not copy and perturb initially, although the
127 // paper does. We follow the code as it should match their reported results
128 // and likely does a better job bouncing out of local optima.
129 // GeneratePerturbedX(X, X);
130
131 // Initialize the unmixing matrix to the whitening matrix.
132 Timer::Start("radical_do_radical");
133 matW = matWhitening;
134
135 mat matYSubspace(nPoints, 2);
136
137 mat matJ = eye(nDims, nDims);
138
139 for (size_t sweepNum = 0; sweepNum < sweeps; sweepNum++)
140 {
141 Log::Info << "RADICAL: sweep " << sweepNum << "." << std::endl;
142
143 for (size_t i = 0; i < nDims - 1; ++i)
144 {
145 for (size_t j = i + 1; j < nDims; ++j)
146 {
147 Log::Debug << "RADICAL 2D on dimensions " << i << " and " << j << "."
148 << std::endl;
149
150 matYSubspace.col(0) = matY.col(i);
151 matYSubspace.col(1) = matY.col(j);
152
153 const double thetaOpt = DoRadical2D(matYSubspace);
154
155 const double cosThetaOpt = cos(thetaOpt);
156 const double sinThetaOpt = sin(thetaOpt);
157
158 // Set elements of transformation matrix.
159 matJ(i, i) = cosThetaOpt;
160 matJ(j, i) = -sinThetaOpt;
161 matJ(i, j) = sinThetaOpt;
162 matJ(j, j) = cosThetaOpt;
163
164 matY *= matJ;
165
166 // Unset elements of transformation matrix, so matJ = eye() again.
167 matJ(i, i) = 1;
168 matJ(j, i) = 0;
169 matJ(i, j) = 0;
170 matJ(j, j) = 1;
171 }
172 }
173 }
174 Timer::Stop("radical_do_radical");
175
176 // The final transposes provide W and Y in the typical form from the ICA
177 // literature.
178 Timer::Start("radical_transpose_data");
179 matW = trans(matW);
180 matY = trans(matY);
181 Timer::Stop("radical_transpose_data");
182 }
183
WhitenFeatureMajorMatrix(const mat & matX,mat & matXWhitened,mat & matWhitening)184 void mlpack::radical::WhitenFeatureMajorMatrix(const mat& matX,
185 mat& matXWhitened,
186 mat& matWhitening)
187 {
188 mat matU, matV;
189 vec s;
190 arma::svd(matU, s, matV, cov(matX));
191 matWhitening = matU * diagmat(1 / sqrt(s)) * trans(matV);
192 matXWhitened = matX * matWhitening;
193 }
194