Update pvl_lambertw.m for faster convergence #18
Merged
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
The precision of the convergence test is too small. As a result you can have vector elements alternately meeting or not meeting the convergence criterion; and since these are all evaluated together, it is common to terminate the loop only after the maximum number of iterations. The attached code demonstrates this.
% pvl_lambertw_test
% Jonathan Allen 2018-12-18
%% set up test
% start with large number of random numbers to find a few that do not
% converge to machine precision
rng('default'); % initialize seed
z = 10 * rand(1e3,1);
% keep track of convergence
w2 = zeros(length(z),36);
converge = true(length(z),36);
%% run pvl_lambertw
% Use asymptotic expansion w = log(z) - log(log(z)) for most z
tmp = log(z + (z == 0));
w = tmp - log(tmp + (tmp == 0));
% Use a series expansion when close to the branch point -1/e
k = (abs(z + 0.3678794411714423216) <= 1.5);
tmp = sqrt(5.43656365691809047*z + 2) - 1; % [1], Eq. 4.22 and text
w(k) = tmp(k);
for k = 1:36
% Converge with Halley's method ([1], Eq. 5.9), about 5 iterations satisfies
% the tolerance for most z
c1 = exp(w);
c2 = w.*c1 - z;
w1 = w + (w ~= -1);
dw = c2./(c1.*w1 - ((w + 2).*c2./(2*w1)));
w = w - dw;
w2(:,k) = w;
converge(:,k) = abs(dw) < 0.7e-16*(2+abs(w));
if all(abs(dw) < 0.7e-16*(2+abs(w)))
%* if all(abs(dw) < eps*(2+abs(w))) %* recommended fix
break;
end
end
%% analyze test
figure
imagesc(converge);
% non convergent values
idx = find(~converge(:,35) | ~converge(:,36));
for i = 1:length(idx)
fprintf('z = %f, w = %f \n',z(idx(i)),w(idx(i)));
end
Collaborator
|
Can you correct this line in your test code?
|
Author
|
Cliff,
I have
dw = c2./(c1.*w1 - ((w + 2).*c2./(2*w1)));
in my test code. Maybe the asterisks got lost in the GitHub comment. It wouldn’t run as you have it. I’ll post as a new function to avoid this.
I’m running on MATLAB 2018b. If you have a different version, you may get different output. Here is mine…
> pvl_lambertw_test
z = 4.018080, w = 1.204631
z = 8.003306, w = 1.606067
z = 8.177606, w = 1.619365
z = 4.161586, w = 1.223875
z = 4.243095, w = 1.234570
z = 8.540999, w = 1.646330
z = 4.053154, w = 1.209384
And you may have to increase the number of random numbers to reproduce the effect.
Thanks,
Jon
From: Cliff Hansen [mailto:notifications@github.com]
Sent: Tuesday, December 18, 2018 1:52 PM
To: sandialabs/MATLAB_PV_LIB <MATLAB_PV_LIB@noreply.github.com>
Cc: Jonathan Allen <jon@allen-analytics.com>; Author <author@noreply.github.com>
Subject: Re: [sandialabs/MATLAB_PV_LIB] Update pvl_lambertw.m for faster convergence (#18)
Can you correct this line in your test code?
dw = c2./(c1.*w1 - ((w + 2).c2./(2w1)));
—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub<#18 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ATCEjZq1nNIcVuYbR0k9lVWhFiVT84tBks5u6VVlgaJpZM4ZZDOh>.
|
cwhanse
reviewed
Dec 18, 2018
Author
|
Yea, I only see it when I pass vectors to pvl_singlediode. The change sped up pvl_singlediode by ~4x for me.
I expect to be back in ABQ in May. Hope to see you then.
From: Cliff Hansen [mailto:notifications@github.com]
Sent: Tuesday, December 18, 2018 2:13 PM
To: sandialabs/MATLAB_PV_LIB <MATLAB_PV_LIB@noreply.github.com>
Cc: Jonathan Allen <jon@allen-analytics.com>; Author <author@noreply.github.com>
Subject: Re: [sandialabs/MATLAB_PV_LIB] Update pvl_lambertw.m for faster convergence (#18)
@cwhanse commented on this pull request.
________________________________
In pvl_lambertw.m<#18 (comment)>:
@@ -55,7 +55,7 @@
dw = c2./(c1.*w1 - ((w + 2).*c2./(2*w1)));
w = w - dw;
- if all(abs(dw) < 0.7e-16*(2+abs(w)))
+ if all(abs(dw) < eps*(2+abs(w)))
Oh nvm, that value illustrates your point. We're cycling in and out of machine precision. Good catch and thanks for the PR!
—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub<#18 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ATCEjXwYw-PHLhAkZiNOvOFZ77YM-TAKks5u6VpqgaJpZM4ZZDOh>.
|
Collaborator
|
thanks very much @joallen6 first PR merged! |
Collaborator
Author
|
Thanks. Just change ‘John’ to ‘Jonathan’.
From: Cliff Hansen [mailto:notifications@github.com]
Sent: Tuesday, December 18, 2018 2:31 PM
To: sandialabs/MATLAB_PV_LIB <MATLAB_PV_LIB@noreply.github.com>
Cc: Jonathan Allen <jon@allen-analytics.com>; Mention <mention@noreply.github.com>
Subject: Re: [sandialabs/MATLAB_PV_LIB] Update pvl_lambertw.m for faster convergence (#18)
@joallen6<https://github.com/joallen6> are you OK with the content of whatsnew<https://github.com/sandialabs/MATLAB_PV_LIB/blob/master/whatsnew.rst>
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub<#18 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ATCEjWYzB63E6rTCD-oJfuSKXwI8zZcFks5u6V6ZgaJpZM4ZZDOh>.
|
Collaborator
|
Done, and thanks again!
From: joallen6 <notifications@github.com>
Sent: Tuesday, December 18, 2018 2:33 PM
To: sandialabs/MATLAB_PV_LIB <MATLAB_PV_LIB@noreply.github.com>
Cc: Hansen, Clifford W <cwhanse@sandia.gov>; State change <state_change@noreply.github.com>
Subject: [EXTERNAL] Re: [sandialabs/MATLAB_PV_LIB] Update pvl_lambertw.m for faster convergence (#18)
Thanks. Just change ‘John’ to ‘Jonathan’.
From: Cliff Hansen [mailto:notifications@github.com]
Sent: Tuesday, December 18, 2018 2:31 PM
To: sandialabs/MATLAB_PV_LIB <MATLAB_PV_LIB@noreply.github.com<mailto:MATLAB_PV_LIB@noreply.github.com>>
Cc: Jonathan Allen <jon@allen-analytics.com<mailto:jon@allen-analytics.com>>; Mention <mention@noreply.github.com<mailto:mention@noreply.github.com>>
Subject: Re: [sandialabs/MATLAB_PV_LIB] Update pvl_lambertw.m for faster convergence (#18)
@joallen6<https://github.com/joallen6> are you OK with the content of whatsnew<https://github.com/sandialabs/MATLAB_PV_LIB/blob/master/whatsnew.rst>
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub<#18 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ATCEjWYzB63E6rTCD-oJfuSKXwI8zZcFks5u6V6ZgaJpZM4ZZDOh>.
—
You are receiving this because you modified the open/close state.
Reply to this email directly, view it on GitHub<#18 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AFJNLzUmftJLD2DRFHf7t_WEBiQvLlnnks5u6V8pgaJpZM4ZZDOh>.
|
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
The precision of the convergence test is too small. As a result you can have vector elements alternately meeting or not meeting the convergence criterion; and since these are all evaluated together, it is common to terminate the loop only after the maximum number of iterations. The attached code demonstrates this.
% pvl_lambertw_test
% Jonathan Allen 2018-12-18
%% set up test
% start with large number of random numbers to find a few that do not
% converge to machine precision
rng('default'); % initialize seed
z = 10 * rand(1e3,1);
% keep track of convergence
w2 = zeros(length(z),36);
converge = true(length(z),36);
%% run pvl_lambertw
% Use asymptotic expansion w = log(z) - log(log(z)) for most z
tmp = log(z + (z == 0));
w = tmp - log(tmp + (tmp == 0));
% Use a series expansion when close to the branch point -1/e
k = (abs(z + 0.3678794411714423216) <= 1.5);
tmp = sqrt(5.43656365691809047*z + 2) - 1; % [1], Eq. 4.22 and text
w(k) = tmp(k);
for k = 1:36
% Converge with Halley's method ([1], Eq. 5.9), about 5 iterations satisfies
% the tolerance for most z
c1 = exp(w);
c2 = w.*c1 - z;
w1 = w + (w ~= -1);
dw = c2./(c1.*w1 - ((w + 2).c2./(2w1)));
w = w - dw;
w2(:,k) = w;
converge(:,k) = abs(dw) < 0.7e-16*(2+abs(w));
if all(abs(dw) < 0.7e-16*(2+abs(w)))
%* if all(abs(dw) < eps*(2+abs(w))) %* recommended fix
break;
end
end
%% analyze test
figure
imagesc(converge);
% non convergent values
idx = find(~converge(:,35) | ~converge(:,36));
for i = 1:length(idx)
fprintf('z = %f, w = %f \n',z(idx(i)),w(idx(i)));
end