fixed performance issue of QR deomposition if #rows >> #columns #54

Merged
merged 1 commit into from Feb 26, 2016
Jump to file or symbol
Failed to load files and symbols.
+59 −19
Diff settings

Always

Just for now

@@ -39,7 +39,6 @@ namespace QuantLib {
MINPACK::qrfac(m, n, mT.begin(), 0, (pivot)?1:0,
lipvt.get(), n, rdiag.get(), rdiag.get(), wa.get());
-
if (r.columns() != n || r.rows() !=n)
r = Matrix(n, n);
@@ -58,27 +57,58 @@ namespace QuantLib {
if (q.rows() != m || q.columns() != n)
q = Matrix(m, n);
- Array w(m);
- for (Size k=0; k < m; ++k) {
- std::fill(w.begin(), w.end(), 0.0);
- w[k] = 1.0;
-
- for (Size j=0; j < std::min(n, m); ++j) {
- const Real t3 = mT[j][j];
- if (t3 != 0.0) {
- const Real t
- = std::inner_product(mT.row_begin(j)+j, mT.row_end(j),
- w.begin()+j, 0.0)/t3;
- for (Size i=j; i<m; ++i) {
- w[i]-=mT[j][i]*t;
+ if (m > n) {
+ std::fill(q.begin(), q.end(), 0.0);
+
+ Integer u = std::min(n,m);
+ for (Size i=0; i < u; ++i)
+ q[i][i] = 1.0;
+
+ Array v(m);
+ for (Integer i=u-1; i >=0; --i) {
+ if (std::fabs(mT[i][i]) > QL_EPSILON) {
+ const Real tau = 1.0/mT[i][i];
+
+ std::fill(v.begin(), v.begin()+i, 0.0);
+ std::copy(mT.row_begin(i)+i, mT.row_end(i), v.begin()+i);
+
+ Array w(n, 0.0);
+ for (Size l=0; l < n; ++l)
+ w[l] += std::inner_product(
+ v.begin()+i, v.end(), q.column_begin(l)+i, 0.0);
+
+ for (Size k=i; k < m; ++k) {
+ const Real a = tau*v[k];
+ for (Size l=0; l < n; ++l)
+ q[k][l] -= a*w[l];
+ }
+ }
+ }
+ }
+ else {
+ Array w(m);
+ for (Size k=0; k < m; ++k) {
+ std::fill(w.begin(), w.end(), 0.0);
+ w[k] = 1.0;
+
+ for (Size j=0; j < std::min(n, m); ++j) {
+ const Real t3 = mT[j][j];
+ if (t3 != 0.0) {
+ const Real t
+ = std::inner_product(mT.row_begin(j)+j, mT.row_end(j),
+ w.begin()+j, 0.0)/t3;
+ for (Size i=j; i<m; ++i) {
+ w[i]-=mT[j][i]*t;
+ }
}
+ q[k][j] = w[j];
}
- q[k][j] = w[j];
+ std::fill(q.row_begin(k) + std::min(n, m), q.row_end(k), 0.0);
}
- std::fill(q.row_begin(k) + std::min(n, m), q.row_end(k), 0.0);
}
std::vector<Size> ipvt(n);
+
if (pivot) {
std::copy(lipvt.get(), lipvt.get()+n, ipvt.begin());
}
@@ -102,10 +132,10 @@ namespace QuantLib {
Matrix q(m, n), r(n, n);
std::vector<Size> lipvt = qrDecomposition(a, q, r, pivot);
+
boost::scoped_array<int> ipvt(new int[n]);
std::copy(lipvt.begin(), lipvt.end(), ipvt.get());
- Matrix aT = transpose(a);
Matrix rT = transpose(r);
boost::scoped_array<Real> sdiag(new Real[n]);
@@ -115,6 +145,7 @@ namespace QuantLib {
if (!d.empty()) {
std::copy(d.begin(), d.end(), ld.begin());
}
+
Array x(n);
Array qtb = transpose(q)*b;
View
@@ -237,6 +237,7 @@ void MatricesTest::testQRDecomposition() {
Matrix Q, R;
bool pivot = true;
const Matrix& A = testMatrices[j];
+
const std::vector<Size> ipvt = qrDecomposition(A, Q, R, pivot);
Matrix P(A.columns(), A.columns(), 0.0);
@@ -271,8 +272,16 @@ void MatricesTest::testQRSolve() {
for (Size i=0; i < std::min(bigM.rows(), bigM.columns()); ++i) {
bigM[i][i] = i+1.0;
}
- Matrix testMatrices[] = { M1, M2, M3, transpose(M3),
- M4, transpose(M4), M5, I, M7, bigM, transpose(bigM) };
+
+ Matrix randM(50, 200);
+ for (Size i=0; i < randM.rows(); ++i)
+ for (Size j=0; j < randM.columns(); ++j)
+ randM[i][j] = rng.next().value;
+
+ Matrix testMatrices[] = {M1, M2, M3, transpose(M3),
+ M4, transpose(M4), M5, I, M7,
+ bigM, transpose(bigM),
+ randM, transpose(randM) };
for (Size j = 0; j < LENGTH(testMatrices); j++) {
const Matrix& A = testMatrices[j];