Skip to content

Fixes for structured LBFGS approximation. #38

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
lambe wants to merge 5 commits into
base: develop
Choose a base branch
from
Open

Fixes for structured LBFGS approximation. #38

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9613f3d
Bugfixes to Structured LBFGS approximation.
Dec 13, 2016
b534954
Refactored structured LBFGS matvec function to fix a major bug.
Dec 15, 2016
91b1d30
Modified structured LBFGS operator to suppress Python runtime
Jan 13, 2017
11d9611
Major bugfix and code cleanup in structured LSR1 operator.
Jan 14, 2017
321f418
Refactored structured LBFGS operator for clarity in the matvec.
Jan 14, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
43 changes: 27 additions & 16 deletions pykrylov/linop/lbfgs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -263,20 +263,24 @@ 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)
self.accept_threshold = 1e-8
super(StructuredLBFGSOperator, self).__init__(n, npairs, **kwargs)
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.matvec(new_s)
ypBs = ys + (ys * np.dot(new_s, Bs))**0.5
Bs = self.qn_matvec(new_s)
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):
Expand Down Expand Up @@ -311,22 +315,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
# Form a[] and ad[] vectors for current step
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:
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])
update = (aTsk / aTsl) * ad[:, l] + (adTsk / aTsl) * a[:, l] - \
(aTsk * adTsl / aTsl**2) * a[:, l]
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])
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
22 changes: 14 additions & 8 deletions pykrylov/linop/lsr1.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Copy link
Collaborator

@dpo dpo Mar 30, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wasn't this already computed in aTs[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()
Copy link
Collaborator

@dpo dpo Mar 30, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why these changes? It's hard to tell whether they actually modify the behavior or are just a rewrite.

Copy link
Author

@lambe lambe Mar 30, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I remember correctly, it's just a rewrite to make the similarity between structured LBFGS and structured LSR1 obvious.

Copy link
Collaborator

@dpo dpo Mar 30, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok. In that case, it's wasteful to recompute the dot products that are already stored in aTs.

Copy link
Author

@lambe lambe Mar 30, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I'll simplify it.


# 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