From 9613f3d0bb2da11fd774898b631e8874b57c1a37 Mon Sep 17 00:00:00 2001 From: Andrew Lambe Date: Tue, 13 Dec 2016 16:05:59 -0500 Subject: [PATCH 1/5] Bugfixes to Structured LBFGS approximation. --- pykrylov/linop/lbfgs.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pykrylov/linop/lbfgs.py b/pykrylov/linop/lbfgs.py index b6f3b3b..3d16bee 100644 --- a/pykrylov/linop/lbfgs.py +++ b/pykrylov/linop/lbfgs.py @@ -263,8 +263,7 @@ def __init__(self, n, npairs=5, **kwargs): Nocedal; the scaling factor is sᵀy/yᵀy (default: False). """ - super(StructuredLBFGSOperator, self).__init__(self, n, npairs, - **kwargs) + super(StructuredLBFGSOperator, self).__init__(n, npairs, **kwargs) self.accept_threshold = 1e-8 def _storing_test(self, new_s, new_y, new_yd, ys): @@ -274,7 +273,7 @@ def _storing_test(self, new_s, new_y, new_yd, ys): ∣yᵀs + √(yᵀs sᵀBs)∣ ⩾ 1e-8. """ - Bs = self.matvec(new_s) + Bs = self.qn_matvec(new_s) ypBs = ys + (ys * np.dot(new_s, Bs))**0.5 return ypBs >= self.accept_threshold From b53495465045534a21d64b96b73d379537bc1364 Mon Sep 17 00:00:00 2001 From: Andrew Lambe Date: Thu, 15 Dec 2016 14:42:58 -0500 Subject: [PATCH 2/5] Refactored structured LBFGS matvec function to fix a major bug. The function now produces the same result as the other LBFGS approximations when y and yd are identical vectors. --- pykrylov/linop/lbfgs.py | 29 ++++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/pykrylov/linop/lbfgs.py b/pykrylov/linop/lbfgs.py index 3d16bee..e047b6c 100644 --- a/pykrylov/linop/lbfgs.py +++ b/pykrylov/linop/lbfgs.py @@ -310,22 +310,29 @@ def qn_matvec(self, v): for i in range(npairs): k = (self.insert + i) % npairs if ys[k] is not None: - coef = (self.gamma * ys[k] / np.dot(s[:, k], s[:, k]))**0.5 - a[:, k] = y[:, k] + coef * s[:, k] / self.gamma - ad[:, k] = yd[:, k] - s[:, k] / self.gamma + # Form a[] and ad[] vectors for current step + a[:, k] = y[:, k].copy() + ad[:, k] = yd[:, k].copy() + Bs = s[:, k] / self.gamma for j in range(i): l = (self.insert + j) % npairs if ys[l] is not None: - alTs = np.dot(a[:, l], s[:, k]) / aTs[l] - adlTs = np.dot(ad[:, l], s[:, k]) - update = alTs / aTs[l] * ad[:, l] + adlTs / aTs[l] * \ - a[:, l] - adTs[l] / aTs[l] * alTs * a[:, l] - a[:, k] += coef * update - ad[:, k] -= update + aTsk = np.dot(a[:, l], s[:, k]) + adTsk = np.dot(ad[:, l], s[:, k]) + aTsl = np.dot(a[:, l], s[:, l]) + adTsl = np.dot(ad[:, l], s[:, l]) + Bs += (aTsk / aTsl) * ad[:, l] + (adTsk / aTsl) * a[:, l] - \ + (aTsk * adTsl / aTsl**2) * a[:, l] + a[:, k] += (ys[k] / np.dot(s[:, k], Bs))**0.5 * Bs + ad[:, k] -= Bs + + # Form inner products with current s[] and input vector aTs[k] = np.dot(a[:, k], s[:, k]) adTs[k] = np.dot(ad[:, k], s[:, k]) aTv = np.dot(a[:, k], v[:]) adTv = np.dot(ad[:, k], v[:]) - q += aTv / aTs[k] * ad[:, k] + adTv / aTs[k] * \ - a[:, k] - aTv * adTs[k] / aTs[k]**2 * a[:, k] + + q += (aTv / aTs[k]) * ad[:, k] + (adTv / aTs[k]) * a[:, k] - \ + (aTv * adTs[k] / aTs[k]**2) * a[:, k] + return q From 91b1d30d2bbd847dbc7fc4f307d3576b96a3bda9 Mon Sep 17 00:00:00 2001 From: Andrew Lambe Date: Fri, 13 Jan 2017 10:44:13 -0500 Subject: [PATCH 3/5] Modified structured LBFGS operator to suppress Python runtime warnings about negative values in the square root and changed the step acceptance threshold to be more permissive. --- pykrylov/linop/lbfgs.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/pykrylov/linop/lbfgs.py b/pykrylov/linop/lbfgs.py index e047b6c..939389a 100644 --- a/pykrylov/linop/lbfgs.py +++ b/pykrylov/linop/lbfgs.py @@ -264,18 +264,23 @@ def __init__(self, n, npairs=5, **kwargs): (default: False). """ super(StructuredLBFGSOperator, self).__init__(n, npairs, **kwargs) - self.accept_threshold = 1e-8 + self.accept_threshold = 1e-10 def _storing_test(self, new_s, new_y, new_yd, ys): u"""Test if new pair {s, y, yd} is to be stored. A new pair {s, y, yd} is only accepted if - ∣yᵀs + √(yᵀs sᵀBs)∣ ⩾ 1e-8. + ∣yᵀs + √(yᵀs sᵀBs)∣ ⩾ self.accept_threshold """ Bs = self.qn_matvec(new_s) - ypBs = ys + (ys * np.dot(new_s, Bs))**0.5 + sBs = np.dot(new_s, Bs) + # Supress python runtime warnings + if ys < 0.0 or sBs < 0.0: + return False + + ypBs = ys + (ys * sBs)**0.5 return ypBs >= self.accept_threshold def qn_matvec(self, v): From 11d9611227310c704c0a5dc3f9093d86f5e28eff Mon Sep 17 00:00:00 2001 From: Andrew Lambe Date: Sat, 14 Jan 2017 11:19:10 -0500 Subject: [PATCH 4/5] Major bugfix and code cleanup in structured LSR1 operator. --- pykrylov/linop/lsr1.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/pykrylov/linop/lsr1.py b/pykrylov/linop/lsr1.py index 0273073..e11a1df 100644 --- a/pykrylov/linop/lsr1.py +++ b/pykrylov/linop/lsr1.py @@ -349,21 +349,27 @@ def qn_matvec(self, v): for i in range(npairs): k = (self.insert + i) % npairs if ys[k] is not None: + # Form all a[] and ad[] vectors for the current step a[:, k] = y[:, k] - s[:, k] / self.gamma ad[:, k] = yd[:, k] - s[:, k] / self.gamma for j in range(i): l = (self.insert + j) % npairs if ys[l] is not None: - alTs = np.dot(a[:, l], s[:, k]) - adlTs = np.dot(ad[:, l], s[:, k]) - update = -alTs / aTs[l] * ad[:, l] - adlTs / aTs[l] * \ - a[:, l] + adTs[l] / aTs[l] * alTs * a[:, l] - a[:, k] += update.copy() - ad[:, k] += update.copy() + aTsk = np.dot(a[:, l], s[:, k]) + adTsk = np.dot(ad[:, l], s[:, k]) + aTsl = np.dot(a[:, l], s[:, l]) + adTsl = np.dot(ad[:, l], s[:, l]) + update = (aTsk / aTsl) * ad[:, l] + (adTsk / aTsl) * a[:, l] - \ + (aTsk * adTsl / aTsl**2) * a[:, l] + a[:, k] -= update.copy() + ad[:, k] -= update.copy() + + # Form inner products with current s[] and input vector aTs[k] = np.dot(a[:, k], s[:, k]) adTs[k] = np.dot(ad[:, k], s[:, k]) aTv = np.dot(a[:, k], v[:]) adTv = np.dot(ad[:, k], v[:]) - q += aTv / aTs[k] * ad[:, k] + adTv / aTs[k] * \ - a[:, k] - aTv * adTs[k] / aTs[k]**2 * a[:, k] + + q += (aTv / aTs[k]) * ad[:, k] + (adTv / aTs[k]) * a[:, k] - \ + (aTv * adTs[k] / aTs[k]**2) * a[:, k] return q From 321f418d37702e6084ec3e169dadb4159006603a Mon Sep 17 00:00:00 2001 From: Andrew Lambe Date: Sat, 14 Jan 2017 13:41:15 -0500 Subject: [PATCH 5/5] Refactored structured LBFGS operator for clarity in the matvec. --- pykrylov/linop/lbfgs.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pykrylov/linop/lbfgs.py b/pykrylov/linop/lbfgs.py index 939389a..ed6d192 100644 --- a/pykrylov/linop/lbfgs.py +++ b/pykrylov/linop/lbfgs.py @@ -316,9 +316,8 @@ def qn_matvec(self, v): k = (self.insert + i) % npairs if ys[k] is not None: # Form a[] and ad[] vectors for current step - a[:, k] = y[:, k].copy() - ad[:, k] = yd[:, k].copy() - Bs = s[:, k] / self.gamma + ad[:, k] = yd[:, k] - s[:, k] / self.gamma + Bsk = s[:, k] / self.gamma for j in range(i): l = (self.insert + j) % npairs if ys[l] is not None: @@ -326,10 +325,11 @@ def qn_matvec(self, v): adTsk = np.dot(ad[:, l], s[:, k]) aTsl = np.dot(a[:, l], s[:, l]) adTsl = np.dot(ad[:, l], s[:, l]) - Bs += (aTsk / aTsl) * ad[:, l] + (adTsk / aTsl) * a[:, l] - \ + update = (aTsk / aTsl) * ad[:, l] + (adTsk / aTsl) * a[:, l] - \ (aTsk * adTsl / aTsl**2) * a[:, l] - a[:, k] += (ys[k] / np.dot(s[:, k], Bs))**0.5 * Bs - ad[:, k] -= Bs + Bsk += update.copy() + ad[:, k] -= update.copy() + a[:, k] = y[:, k] + (ys[k] / np.dot(s[:, k], Bsk))**0.5 * Bsk # Form inner products with current s[] and input vector aTs[k] = np.dot(a[:, k], s[:, k])