Skip to content

O
opencv
Commit 5368a4ac authored by Alexander Alekhin's avatar Alexander Alekhin
Browse files
Options

Merge pull request #14102 from alalek:core_refactor_eigenvalues

parents 9340fc0c 5451b89a
Show whitespace changes
Inline Side-by-side
Showing with 102 additions and 95 deletions
modules/core/src/lda.cpp
View file @ 5368a4ac
...@@ -248,9 +248,6 @@ private: ...@@ -248,9 +248,6 @@ private:
// Holds the data dimension. // Holds the data dimension.
int n; int n;
// Stores real/imag part of a complex division.
double cdivr, cdivi;
// Pointer to internal memory. // Pointer to internal memory.
double *d, *e, *ort; double *d, *e, *ort;
double **V, **H; double **V, **H;
...@@ -297,8 +294,9 @@ private: ...@@ -297,8 +294,9 @@ private:
return arr; return arr;
} }
void cdiv(double xr, double xi, double yr, double yi) { static void complex_div(double xr, double xi, double yr, double yi, double& cdivr, double& cdivi) {
double r, dv; double r, dv;
CV_DbgAssert(std::abs(yr) + std::abs(yi) > 0.0);
if (std::abs(yr) > std::abs(yi)) { if (std::abs(yr) > std::abs(yi)) {
r = yi / yr; r = yi / yr;
dv = yr + r * yi; dv = yr + r * yi;
...@@ -324,24 +322,25 @@ private: ...@@ -324,24 +322,25 @@ private:
// Initialize // Initialize
const int max_iters_count = 1000 * this->n; const int max_iters_count = 1000 * this->n;
int nn = this->n; const int nn = this->n; CV_Assert(nn > 0);
int n1 = nn - 1; int n1 = nn - 1;
int low = 0; const int low = 0;
int high = nn - 1; const int high = nn - 1;
double eps = std::pow(2.0, -52.0); const double eps = std::numeric_limits<double>::epsilon();
double exshift = 0.0; double exshift = 0.0;
double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;
// Store roots isolated by balanc and compute matrix norm // Store roots isolated by balanc and compute matrix norm
double norm = 0.0; double norm = 0.0;
for (int i = 0; i < nn; i++) { for (int i = 0; i < nn; i++) {
#if 0 // 'if' condition is always false
if (i < low || i > high) { if (i < low || i > high) {
d[i] = H[i][i]; d[i] = H[i][i];
e[i] = 0.0; e[i] = 0.0;
} }
#endif
for (int j = std::max(i - 1, 0); j < nn; j++) { for (int j = std::max(i - 1, 0); j < nn; j++) {
norm = norm + std::abs(H[i][j]); norm += std::abs(H[i][j]);
} }
} }
...@@ -355,7 +354,7 @@ private: ...@@ -355,7 +354,7 @@ private:
if (norm < FLT_EPSILON) { if (norm < FLT_EPSILON) {
break; break;
} }
s = std::abs(H[l - 1][l - 1]) + std::abs(H[l][l]); double s = std::abs(H[l - 1][l - 1]) + std::abs(H[l][l]);
if (s == 0.0) { if (s == 0.0) {
s = norm; s = norm;
} }
...@@ -366,29 +365,26 @@ private: ...@@ -366,29 +365,26 @@ private:
} }
// Check for convergence // Check for convergence
// One root found
if (l == n1) { if (l == n1) {
// One root found
H[n1][n1] = H[n1][n1] + exshift; H[n1][n1] = H[n1][n1] + exshift;
d[n1] = H[n1][n1]; d[n1] = H[n1][n1];
e[n1] = 0.0; e[n1] = 0.0;
n1--; n1--;
iter = 0; iter = 0;
// Two roots found
} else if (l == n1 - 1) { } else if (l == n1 - 1) {
w = H[n1][n1 - 1] * H[n1 - 1][n1]; // Two roots found
p = (H[n1 - 1][n1 - 1] - H[n1][n1]) / 2.0; double w = H[n1][n1 - 1] * H[n1 - 1][n1];
q = p * p + w; double p = (H[n1 - 1][n1 - 1] - H[n1][n1]) * 0.5;
z = std::sqrt(std::abs(q)); double q = p * p + w;
double z = std::sqrt(std::abs(q));
H[n1][n1] = H[n1][n1] + exshift; H[n1][n1] = H[n1][n1] + exshift;
H[n1 - 1][n1 - 1] = H[n1 - 1][n1 - 1] + exshift; H[n1 - 1][n1 - 1] = H[n1 - 1][n1 - 1] + exshift;
x = H[n1][n1]; double x = H[n1][n1];
// Real pair
if (q >= 0) { if (q >= 0) {
// Real pair
if (p >= 0) { if (p >= 0) {
z = p + z; z = p + z;
} else { } else {
...@@ -402,10 +398,10 @@ private: ...@@ -402,10 +398,10 @@ private:
e[n1 - 1] = 0.0; e[n1 - 1] = 0.0;
e[n1] = 0.0; e[n1] = 0.0;
x = H[n1][n1 - 1]; x = H[n1][n1 - 1];
s = std::abs(x) + std::abs(z); double s = std::abs(x) + std::abs(z);
p = x / s; p = x / s;
q = z / s; q = z / s;
r = std::sqrt(p * p + q * q); double r = std::sqrt(p * p + q * q);
p = p / r; p = p / r;
q = q / r; q = q / r;
...@@ -433,9 +429,8 @@ private: ...@@ -433,9 +429,8 @@ private:
V[i][n1] = q * V[i][n1] - p * z; V[i][n1] = q * V[i][n1] - p * z;
} }
// Complex pair
} else { } else {
// Complex pair
d[n1 - 1] = x + p; d[n1 - 1] = x + p;
d[n1] = x + p; d[n1] = x + p;
e[n1 - 1] = z; e[n1 - 1] = z;
...@@ -444,28 +439,25 @@ private: ...@@ -444,28 +439,25 @@ private:
n1 = n1 - 2; n1 = n1 - 2;
iter = 0; iter = 0;
// No convergence yet
} else { } else {
// No convergence yet
// Form shift // Form shift
double x = H[n1][n1];
x = H[n1][n1]; double y = 0.0;
y = 0.0; double w = 0.0;
w = 0.0;
if (l < n1) { if (l < n1) {
y = H[n1 - 1][n1 - 1]; y = H[n1 - 1][n1 - 1];
w = H[n1][n1 - 1] * H[n1 - 1][n1]; w = H[n1][n1 - 1] * H[n1 - 1][n1];
} }
// Wilkinson's original ad hoc shift // Wilkinson's original ad hoc shift
if (iter == 10) { if (iter == 10) {
exshift += x; exshift += x;
for (int i = low; i <= n1; i++) { for (int i = low; i <= n1; i++) {
H[i][i] -= x; H[i][i] -= x;
} }
s = std::abs(H[n1][n1 - 1]) + std::abs(H[n1 - 1][n1 - 2]); double s = std::abs(H[n1][n1 - 1]) + std::abs(H[n1 - 1][n1 - 2]);
x = y = 0.75 * s; x = y = 0.75 * s;
w = -0.4375 * s * s; w = -0.4375 * s * s;
} }
...@@ -473,14 +465,14 @@ private: ...@@ -473,14 +465,14 @@ private:
// MATLAB's new ad hoc shift // MATLAB's new ad hoc shift
if (iter == 30) { if (iter == 30) {
s = (y - x) / 2.0; double s = (y - x) * 0.5;
s = s * s + w; s = s * s + w;
if (s > 0) { if (s > 0) {
s = std::sqrt(s); s = std::sqrt(s);
if (y < x) { if (y < x) {
s = -s; s = -s;
} }
s = x - w / ((y - x) / 2.0 + s); s = x - w / ((y - x) * 0.5 + s);
for (int i = low; i <= n1; i++) { for (int i = low; i <= n1; i++) {
H[i][i] -= s; H[i][i] -= s;
} }
...@@ -493,12 +485,16 @@ private: ...@@ -493,12 +485,16 @@ private:
if (iter > max_iters_count) if (iter > max_iters_count)
CV_Error(Error::StsNoConv, "Algorithm doesn't converge (complex eigen values?)"); CV_Error(Error::StsNoConv, "Algorithm doesn't converge (complex eigen values?)");
double p = std::numeric_limits<double>::quiet_NaN();
double q = std::numeric_limits<double>::quiet_NaN();
double r = std::numeric_limits<double>::quiet_NaN();
// Look for two consecutive small sub-diagonal elements // Look for two consecutive small sub-diagonal elements
int m = n1 - 2; int m = n1 - 2;
while (m >= l) { while (m >= l) {
z = H[m][m]; double z = H[m][m];
r = x - z; r = x - z;
s = y - z; double s = y - z;
p = (r * s - w) / H[m + 1][m] + H[m][m + 1]; p = (r * s - w) / H[m + 1][m] + H[m][m + 1];
q = H[m + 1][m + 1] - z - r - s; q = H[m + 1][m + 1] - z - r - s;
r = H[m + 2][m + 1]; r = H[m + 2][m + 1];
...@@ -527,6 +523,7 @@ private: ...@@ -527,6 +523,7 @@ private:
// Double QR step involving rows l:n and columns m:n // Double QR step involving rows l:n and columns m:n
for (int k = m; k < n1; k++) { for (int k = m; k < n1; k++) {
bool notlast = (k != n1 - 1); bool notlast = (k != n1 - 1);
if (k != m) { if (k != m) {
p = H[k][k - 1]; p = H[k][k - 1];
...@@ -542,7 +539,7 @@ private: ...@@ -542,7 +539,7 @@ private:
if (x == 0.0) { if (x == 0.0) {
break; break;
} }
s = std::sqrt(p * p + q * q + r * r); double s = std::sqrt(p * p + q * q + r * r);
if (p < 0) { if (p < 0) {
s = -s; s = -s;
} }
...@@ -555,7 +552,7 @@ private: ...@@ -555,7 +552,7 @@ private:
p = p + s; p = p + s;
x = p / s; x = p / s;
y = q / s; y = q / s;
z = r / s; double z = r / s;
q = q / p; q = q / p;
r = r / p; r = r / p;
...@@ -567,8 +564,8 @@ private: ...@@ -567,8 +564,8 @@ private:
p = p + r * H[k + 2][j]; p = p + r * H[k + 2][j];
H[k + 2][j] = H[k + 2][j] - p * z; H[k + 2][j] = H[k + 2][j] - p * z;
} }
H[k][j] = H[k][j] - p * x; H[k][j] -= p * x;
H[k + 1][j] = H[k + 1][j] - p * y; H[k + 1][j] -= p * y;
} }
// Column modification // Column modification
...@@ -579,8 +576,8 @@ private: ...@@ -579,8 +576,8 @@ private:
p = p + z * H[i][k + 2]; p = p + z * H[i][k + 2];
H[i][k + 2] = H[i][k + 2] - p * r; H[i][k + 2] = H[i][k + 2] - p * r;
} }
H[i][k] = H[i][k] - p; H[i][k] -= p;
H[i][k + 1] = H[i][k + 1] - p * q; H[i][k + 1] -= p * q;
} }
// Accumulate transformations // Accumulate transformations
...@@ -606,17 +603,19 @@ private: ...@@ -606,17 +603,19 @@ private:
} }
for (n1 = nn - 1; n1 >= 0; n1--) { for (n1 = nn - 1; n1 >= 0; n1--) {
p = d[n1]; double p = d[n1];
q = e[n1]; double q = e[n1];
if (q == 0) {
// Real vector // Real vector
double z = std::numeric_limits<double>::quiet_NaN();
double s = std::numeric_limits<double>::quiet_NaN();
if (q == 0) {
int l = n1; int l = n1;
H[n1][n1] = 1.0; H[n1][n1] = 1.0;
for (int i = n1 - 1; i >= 0; i--) { for (int i = n1 - 1; i >= 0; i--) {
w = H[i][i] - p; double w = H[i][i] - p;
r = 0.0; double r = 0.0;
for (int j = l; j <= n1; j++) { for (int j = l; j <= n1; j++) {
r = r + H[i][j] * H[j][n1]; r = r + H[i][j] * H[j][n1];
} }
...@@ -631,34 +630,38 @@ private: ...@@ -631,34 +630,38 @@ private:
} else { } else {
H[i][n1] = -r / (eps * norm); H[i][n1] = -r / (eps * norm);
} }
// Solve real equations
} else { } else {
x = H[i][i + 1]; // Solve real equations
y = H[i + 1][i]; CV_DbgAssert(!cvIsNaN(z));
double x = H[i][i + 1];
double y = H[i + 1][i];
q = (d[i] - p) * (d[i] - p) + e[i] * e[i]; q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
t = (x * s - z * r) / q; double t = (x * s - z * r) / q;
H[i][n1] = t; H[i][n1] = t;
if (std::abs(x) > std::abs(z)) { if (std::abs(x) > std::abs(z)) {
H[i + 1][n1] = (-r - w * t) / x; H[i + 1][n1] = (-r - w * t) / x;
} else { } else {
CV_DbgAssert(z != 0.0);
H[i + 1][n1] = (-s - y * t) / z; H[i + 1][n1] = (-s - y * t) / z;
} }
} }
// Overflow control // Overflow control
double t = std::abs(H[i][n1]);
t = std::abs(H[i][n1]);
if ((eps * t) * t > 1) { if ((eps * t) * t > 1) {
double inv_t = 1.0 / t;
for (int j = i; j <= n1; j++) { for (int j = i; j <= n1; j++) {
H[j][n1] = H[j][n1] / t; H[j][n1] *= inv_t;
} }
} }
} }
} }
// Complex vector
} else if (q < 0) { } else if (q < 0) {
// Complex vector
double z = std::numeric_limits<double>::quiet_NaN();
double r = std::numeric_limits<double>::quiet_NaN();
double s = std::numeric_limits<double>::quiet_NaN();
int l = n1 - 1; int l = n1 - 1;
// Last vector component imaginary so matrix is triangular // Last vector component imaginary so matrix is triangular
...@@ -667,9 +670,11 @@ private: ...@@ -667,9 +670,11 @@ private:
H[n1 - 1][n1 - 1] = q / H[n1][n1 - 1]; H[n1 - 1][n1 - 1] = q / H[n1][n1 - 1];
H[n1 - 1][n1] = -(H[n1][n1] - p) / H[n1][n1 - 1]; H[n1 - 1][n1] = -(H[n1][n1] - p) / H[n1][n1 - 1];
} else { } else {
cdiv(0.0, -H[n1 - 1][n1], H[n1 - 1][n1 - 1] - p, q); complex_div(
H[n1 - 1][n1 - 1] = cdivr; 0.0, -H[n1 - 1][n1],
H[n1 - 1][n1] = cdivi; H[n1 - 1][n1 - 1] - p, q,
H[n1 - 1][n1 - 1], H[n1 - 1][n1]
);
} }
H[n1][n1 - 1] = 0.0; H[n1][n1 - 1] = 0.0;
H[n1][n1] = 1.0; H[n1][n1] = 1.0;
...@@ -681,7 +686,7 @@ private: ...@@ -681,7 +686,7 @@ private:
ra = ra + H[i][j] * H[j][n1 - 1]; ra = ra + H[i][j] * H[j][n1 - 1];
sa = sa + H[i][j] * H[j][n1]; sa = sa + H[i][j] * H[j][n1];
} }
w = H[i][i] - p; double w = H[i][i] - p;
if (e[i] < 0.0) { if (e[i] < 0.0) {
z = w; z = w;
...@@ -690,41 +695,42 @@ private: ...@@ -690,41 +695,42 @@ private:
} else { } else {
l = i; l = i;
if (e[i] == 0) { if (e[i] == 0) {
cdiv(-ra, -sa, w, q); complex_div(
H[i][n1 - 1] = cdivr; -ra, -sa,
H[i][n1] = cdivi; w, q,
H[i][n1 - 1], H[i][n1]
);
} else { } else {
// Solve complex equations // Solve complex equations
x = H[i][i + 1]; double x = H[i][i + 1];
y = H[i + 1][i]; double y = H[i + 1][i];
vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q; vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
vi = (d[i] - p) * 2.0 * q; vi = (d[i] - p) * 2.0 * q;
if (vr == 0.0 && vi == 0.0) { if (vr == 0.0 && vi == 0.0) {
vr = eps * norm * (std::abs(w) + std::abs(q) + std::abs(x) vr = eps * norm * (std::abs(w) + std::abs(q) + std::abs(x)
+ std::abs(y) + std::abs(z)); + std::abs(y) + std::abs(z));
} }
cdiv(x * r - z * ra + q * sa, complex_div(
x * s - z * sa - q * ra, vr, vi); x * r - z * ra + q * sa, x * s - z * sa - q * ra,
H[i][n1 - 1] = cdivr; vr, vi,
H[i][n1] = cdivi; H[i][n1 - 1], H[i][n1]);
if (std::abs(x) > (std::abs(z) + std::abs(q))) { if (std::abs(x) > (std::abs(z) + std::abs(q))) {
H[i + 1][n1 - 1] = (-ra - w * H[i][n1 - 1] + q H[i + 1][n1 - 1] = (-ra - w * H[i][n1 - 1] + q
* H[i][n1]) / x; * H[i][n1]) / x;
H[i + 1][n1] = (-sa - w * H[i][n1] - q * H[i][n1 H[i + 1][n1] = (-sa - w * H[i][n1] - q * H[i][n1
- 1]) / x; - 1]) / x;
} else { } else {
cdiv(-r - y * H[i][n1 - 1], -s - y * H[i][n1], z, complex_div(
q); -r - y * H[i][n1 - 1], -s - y * H[i][n1],
H[i + 1][n1 - 1] = cdivr; z, q,
H[i + 1][n1] = cdivi; H[i + 1][n1 - 1], H[i + 1][n1]);
} }
} }
// Overflow control // Overflow control
t = std::max(std::abs(H[i][n1 - 1]), std::abs(H[i][n1])); double t = std::max(std::abs(H[i][n1 - 1]), std::abs(H[i][n1]));
if ((eps * t) * t > 1) { if ((eps * t) * t > 1) {
for (int j = i; j <= n1; j++) { for (int j = i; j <= n1; j++) {
H[j][n1 - 1] = H[j][n1 - 1] / t; H[j][n1 - 1] = H[j][n1 - 1] / t;
...@@ -738,6 +744,7 @@ private: ...@@ -738,6 +744,7 @@ private:
// Vectors of isolated roots // Vectors of isolated roots
#if 0 // 'if' condition is always false
for (int i = 0; i < nn; i++) { for (int i = 0; i < nn; i++) {
if (i < low || i > high) { if (i < low || i > high) {
for (int j = i; j < nn; j++) { for (int j = i; j < nn; j++) {
...@@ -745,14 +752,15 @@ private: ...@@ -745,14 +752,15 @@ private:
} }
} }
} }
#endif
// Back transformation to get eigenvectors of original matrix // Back transformation to get eigenvectors of original matrix
for (int j = nn - 1; j >= low; j--) { for (int j = nn - 1; j >= low; j--) {
for (int i = low; i <= high; i++) { for (int i = low; i <= high; i++) {
z = 0.0; double z = 0.0;
for (int k = low; k <= std::min(j, high); k++) { for (int k = low; k <= std::min(j, high); k++) {
z = z + V[i][k] * H[k][j]; z += V[i][k] * H[k][j];
} }
V[i][j] = z; V[i][j] = z;
} }
...@@ -852,15 +860,15 @@ private: ...@@ -852,15 +860,15 @@ private:
// Releases all internal working memory. // Releases all internal working memory.
void release() { void release() {
// releases the working data // releases the working data
delete[] d; delete[] d; d = NULL;
delete[] e; delete[] e; e = NULL;
delete[] ort; delete[] ort; ort = NULL;
for (int i = 0; i < n; i++) { for (int i = 0; i < n; i++) {
delete[] H[i]; if (H) delete[] H[i];
delete[] V[i]; if (V) delete[] V[i];
} }
delete[] H; delete[] H; H = NULL;
delete[] V; delete[] V; V = NULL;
} }
// Computes the Eigenvalue Decomposition for a matrix given in H. // Computes the Eigenvalue Decomposition for a matrix given in H.
...@@ -870,7 +878,7 @@ private: ...@@ -870,7 +878,7 @@ private:
d = alloc_1d<double> (n); d = alloc_1d<double> (n);
e = alloc_1d<double> (n); e = alloc_1d<double> (n);
ort = alloc_1d<double> (n); ort = alloc_1d<double> (n);
try { {
// Reduce to Hessenberg form. // Reduce to Hessenberg form.
orthes(); orthes();
// Reduce Hessenberg to real Schur form. // Reduce Hessenberg to real Schur form.
...@@ -888,11 +896,6 @@ private: ...@@ -888,11 +896,6 @@ private:
// Deallocate the memory by releasing all internal working data. // Deallocate the memory by releasing all internal working data.
release(); release();
} }
catch (...)
{
release();
throw;
}
} }
public: public:
...@@ -900,7 +903,11 @@ public: ...@@ -900,7 +903,11 @@ public:
// given in src. This function is a port of the EigenvalueSolver in JAMA, // given in src. This function is a port of the EigenvalueSolver in JAMA,
// which has been released to public domain by The MathWorks and the // which has been released to public domain by The MathWorks and the
// National Institute of Standards and Technology (NIST). // National Institute of Standards and Technology (NIST).
EigenvalueDecomposition(InputArray src, bool fallbackSymmetric = true) { EigenvalueDecomposition(InputArray src, bool fallbackSymmetric = true) :
n(0),
d(NULL), e(NULL), ort(NULL),
V(NULL), H(NULL)
{
compute(src, fallbackSymmetric); compute(src, fallbackSymmetric);
} }
...@@ -938,7 +945,7 @@ public: ...@@ -938,7 +945,7 @@ public:
} }
} }
~EigenvalueDecomposition() {} ~EigenvalueDecomposition() { release(); }
// Returns the eigenvalues of the Eigenvalue Decomposition. // Returns the eigenvalues of the Eigenvalue Decomposition.
Mat eigenvalues() const { return _eigenvalues; } Mat eigenvalues() const { return _eigenvalues; }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment