Matlab: Plot objective and gradients (partial derivatives)

Posted by rolf on Thursday, 6 November 2014

Given a function and the gradient/derivative function of it, this code can plot the function and for each dimension the partial derivatives.

If you want to minimize your objective function with a gradient based optimizer, then it sometimes comes handy to visualize the objective functions (or at least projections of it) and the gradient function. Then it is more easy to see if your gradients are correct and if the objective function is optimizable.

Of course another way to check if your gradients are correct is to compute the numerical gradient and compare the difference to your computed gradient.

NOTE ON THE CODE: If your matlab version does not provide the option MaxHeadSize whith quiver, you can have a look at http://www.mathworks.de/matlabcentral/fileexchange/34305-adjust--quiver-....

Here are some example images how the figure look like. In blue the objective function, in green the gradients/tangents.

The CODE for plotting the objective and the gradients:

plotGradAndObj.m

Code (expand / collapse)

function [] = plotGradAndObj(f,X,plotOptions,varargin)
% rolfk, 06.11.2014

% plotGradAndObj plots the objective function around a given startpoint.
% the area [-sStep.*nStep : sStep : sStep.*nStep] around X is plotted.
% Each dimension of X is plotted individually = a projection of each dimension is plotted.
% Dimension m is plotted by varying only x_m and keeping all other x_i fixed
%   to the initialization value.
%
% INPUT:
% % f : a matlab .m file that returns the objective value and the gradient
%     gradient vector of size Nx1     
% X = (x_1,...,x_N)': the point (of size Nx1) where to start plotting

% ---------- plotOptions ----------
% plotOptions.sStep       = scalar or Nx1 vector
%                           determines the distance between points, which are evaluated
%                           on the objective and the gradient
%                           if sStep is a scalar the distance is the same for all variables
%
% plotOptions.nStep       = scalar or Nx1 vector (how many points to evaluate in each variable)
%                           if nStep is a scalar in each direction the same number of points are evaluated
%
% plotOptions.gradStep    = plot gradient arrow every gradArrowStep-th point
%
% plotOptions.dim         = array of size (Mx1) with 1<=M<=N: which dimension to plot
%
% plotOptions.ShowGrad    = 1 (=show Gradient=tangents), 0= do not show gradients
%
% plotOptions.arrowLength = length of tangents in the plot
%
% plotOptions.MaxHeadSize = size of arrow head of the tangents
%
% plotOptions.Marker      = Marker style of the point where the tangent touches the curve

% varargin: if the gradient and objective function require more input it can be provided here.

% REQUIREMENT ON THE FUNCTION f
% [obj, grad] = f(X,aux)
% with obj being the function value and grad the array of partial derivatives.
%
% USAGE - EXAMPLE:
%  plotGradAndObj('exampleObjAndGrad',x0,options,[]);
%
% in a m-file exampleObjAndGrad.m write following lines
% ---------- exampleObjAndGrad.m ----------
%  function [obj, grad] = exampleObjAndGrad(x)
%  % ---------- compute objective ----------
%  x1 = x(1);
%  x2 = x(2);
%  x3 = x(3);
%  x4 = x(4);
%  
%  obj = 1/4*x4^4*x3 + 1/3*x3^3*x2 + 1/2*x2^2 + x1*x3^2 + 1;
%  
%  
%  % ---------- compute gradient ----------
%  grad(1) = x3^2;
%  grad(2) = 1/3*x3^3 + x2;
%  grad(3) = 1/4*x4^4 + x3^2*x2+x1*2*x3;
%  grad(4) = x4^3*x3;
%  
%  return
%
% in another m-file runExample.m following lines:
% ---------- runExample.m ----------
%  x0  = [-1  -1  -1  -1];   % initial point
%  options.sStep = 1e-2;
%  options.nStep = 100;
%  options.gradStep = [];
%  options.dim = [1,2,3,4];
%  
%  % show gradients=tangents or not
%  options.ShowGrad = 1;
%  
%  % how the arrows of the gradients=tangents look
%  options.arrowLength = 0.1;
%  options.MaxHeadSize = 1e-10;
%  options.Marker = 'o';
%  
%  plotGradAndObj('exampleObjAndGrad',x0,options,[]);
% ====================================================================================================

% ---------- PLOTTING PARAMETERS ----------
% optional fields
% see http://www.mathworks.de/help/matlab/ref/quiverseries-properties.html
if (~isfield(plotOptions,'arrowLength') || isempty(plotOptions.arrowLength))
    arrowLength = 0.1;
else
    arrowLength = plotOptions.arrowLength;
end

if (~isfield(plotOptions,'MaxHeadSize') || isempty(plotOptions.MaxHeadSize))
    MaxHeadSize = 1e-10;
else
    MaxHeadSize = plotOptions.MaxHeadSize;
end

if (~isfield(plotOptions,'Marker') || isempty(plotOptions.Marker))
    Marker = 'o'; % 'none'
else
    Marker = plotOptions.Marker;
end

if (~isfield(plotOptions,'sStep') || isempty(plotOptions.sStep))
  sStep = 1e-2;
else
    sStep = plotOptions.sStep;
end

if (~isfield(plotOptions,'nStep') || isempty(plotOptions.nStep))
  nStep = 100;
else
    nStep = plotOptions.nStep;
end

if (~isfield(plotOptions,'gradStep') || isempty(plotOptions.gradStep))
  gradStep = 10;
else
    gradStep = plotOptions.gradStep;
end

if (~isfield(plotOptions,'showGrad') || isempty(plotOptions.showGrad))
  showGrad = 1;
else
    showGrad = plotOptions.showGrad;
end

% required field
if (~isfield(plotOptions,'dim') || isempty(plotOptions.dim))
    error('Please specify the dimensions you want to plot')
else
    dim = plotOptions.dim;
end

if length(nStep) == 1
    nStep = nStep * ones(length(X),1);
end

if length(sStep) == 1
    sStep = sStep * ones(length(X),1);
end

%% compute the point-values for the objective and gradient function
Xorg = X;
for j = 1 : length(dim) % i-th component of vector X
    funVec  = zeros(2*nStep(j)+1,1); % function values at different points
    gradMtx = zeros(2*nStep(j)+1,length(X)); % partial derivatives at different points
    % fprintf('%i of %i \n', i, length(X));
    i = dim(j);
    X = Xorg;
    iiStep = 0;
    X(i) = X(i) - nStep(i)*sStep(i);
    for iStep = -nStep(i): 1: nStep(i)
        iiStep = iiStep + 1;
        X(i) = X(i) + sStep(i);
        if isempty(varargin{:})
            [obj, grad] = feval(f,X);
        else
            [obj, grad] = feval(f,X,varargin{:});
        end
        
        funVec(iiStep)    = obj;
        gradMtx(iiStep,:) = grad;
    end
    
    %% plot the objective and gradient function
    x = -sStep(i)*nStep(i) : sStep(i) : sStep(i)*nStep(i);
    figure, plot(x,funVec,'LineWidth', 2), title(sprintf('Dimension=%i',i));
    
    if showGrad
        hold on;
        xPos = 1 : gradStep : 2*nStep(i)+1;
        xQ = x(xPos); yQ=funVec(xPos)'; vQ=gradMtx(xPos,i)'; uQ=ones(1,length(xPos));
        h = quiver(xQ,yQ,uQ,vQ,arrowLength,'LineWidth',2, 'ShowArrowHead', 'on','MaxHeadSize', MaxHeadSize, 'Marker', Marker );
        hold off;
    end
    drawnow;
end

return
% ----------------------------------------------------------------------------------------------------

% I guess following script is only needed in older matlab versions, where the 'MaxHeadSize' option was not available for quiver
% ====================================================================================================
function adjust_quiver_arrowhead_size(quivergroup_handle, scaling_factor)
% downloaded from http://www.mathworks.de/matlabcentral/fileexchange/34305-adjust--quiver--arrowhead-size
% downloaded at 06.Nov.2014 by rolfk
% accompanies the code plotGradAndObj.m
%
% Copyright (c) 2011, Kevin J. Delaney
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of the BMT Scientific Marine Services nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% Make quiver arrowheads bigger or smaller.
%
% adjust_quiver_arrowhead_size(quivergroup_handle, scaling_factor)
%
% Example:
%   h = quiver(1:100, 1:100, randn(100, 100), randn(100, 100));
%   adjust_quiver_arrowhead_size(h, 1.5);   % Makes all arrowheads 50% bigger.
%
% Inputs:
%   quivergroup_handle      Handle returned by "quiver" command.
%   scaling_factor          Factor by which to shrink/grow arrowheads.
%
% Output: none

% Kevin J. Delaney
% December 21, 2011
% BMT Scientific Marine Services (www.scimar.com)

if ~exist('quivergroup_handle', 'var')
    help(mfilename);
    return
end

if isempty(quivergroup_handle) || any(~ishandle(quivergroup_handle))
    errordlg('Input "quivergroup_handle" is empty or contains invalid handles.', ...
             mfilename);
    return
end

if length(quivergroup_handle) > 1
    errordlg('Expected "quivergroup_handle" to be a single handle.', mfilename);
    return
end

if ~strcmpi(get(quivergroup_handle, 'Type'), 'hggroup')
    errrodlg('Input "quivergroup_handle" is not of type "hggroup".', mfilename);
    return
end

if ~exist('scaling_factor', 'var') || ...
   isempty(scaling_factor) || ...
   ~isnumeric(scaling_factor)
    errordlg('Input "scaling_factor" is missing, empty or non-numeric.', ...
             mfilename);
    return
end

if length(scaling_factor) > 1
    errordlg('Expected "scaling_factor" to be a scalar.', mfilename);
    return
end

if scaling_factor <= 0
    errordlg('"Scaling_factor" should be > 0.', mfilename);
    return
end

line_handles = get(quivergroup_handle, 'Children');

if isempty(line_handles) || (length(line_handles) < 3) || ...
   ~ishandle(line_handles(2)) || ~strcmpi(get(line_handles(2), 'Type'), 'line')
    errordlg('Unable to adjust arrowheads.', mfilename);
    return
end

arrowhead_line = line_handles(2);

XData = get(arrowhead_line, 'XData');
YData = get(arrowhead_line, 'YData');

if isempty(XData) || isempty(YData)
    return
end

%   Break up XData, YData into triplets separated by NaNs.
first_nan_index = find(~isnan(XData), 1, 'first');
last_nan_index  = find(~isnan(XData), 1, 'last');

for index = first_nan_index : 4 : last_nan_index
    these_indices = index + (0:2);
    
    if these_indices(end) > length(XData)
        break
    end
    
    x_triplet = XData(these_indices);
    y_triplet = YData(these_indices);
    
    if any(isnan(x_triplet)) || any(isnan(y_triplet))
        continue
    end
    
    %   First pair.
    delta_x = diff(x_triplet(1:2));
    delta_y = diff(y_triplet(1:2));
    x_triplet(1) = x_triplet(2) - (delta_x * scaling_factor);
    y_triplet(1) = y_triplet(2) - (delta_y * scaling_factor);
        
    %   Second pair.
    delta_x = diff(x_triplet(2:3));
    delta_y = diff(y_triplet(2:3));
    x_triplet(3) = x_triplet(2) + (delta_x * scaling_factor);
    y_triplet(3) = y_triplet(2) + (delta_y * scaling_factor);
    
    XData(these_indices) = x_triplet;
    YData(these_indices) = y_triplet;
end

set(arrowhead_line, 'XData', XData, 'YData', YData);

return
% ----------------------------------------------------------------------------------------------------

And here an example script to call it

RunExample.m

Code (expand / collapse)

% ---------- RunExample.m ----------
x0  = [-1  -1  -1  -1];   % initial point
options.sStep = 1e-2;
options.nStep = 100;
options.gradStep = 5;
options.dim = [1,2,3,4];

% show gradients=tangents or not
options.ShowGrad = 1;

% how the arrows of the gradients=tangents look
options.arrowLength = 0.3;
options.MaxHeadSize = 1e-10;
options.Marker = 'o';

plotGradAndObj('exampleObjAndGrad',x0,options,[]);

And here the Function that returns the objective and gradient:

exampleObjAndGrad.m

Code (expand / collapse)

% ---------- exampleObjAndGrad.m ---------
function [obj, grad] = exampleObjAndGrad(x)
% ---------- compute objective ----------
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);

obj = 1/4*x4^4*x3 + 1/3*x3^3*x2 + 1/2*x2^2 + x1*x3^2 + 1;

% ---------- compute gradient ----------
grad(1) = x3^2;
grad(2) = 1/3*x3^3 + x2;
grad(3) = 1/4*x4^4 + x3^2*x2+x1*2*x3;
grad(4) = x4^3*x3;

return

Files:

blog_tags:

optimization

Add new comment

Your name

Subject

Comment *

More information about text formats

Text format

Filtered HTML

Allowed HTML tags: <a> <cite> <blockquote> <code> <ul> <ol> <li> <dl> <dt> <dd> <pre> <img>
Lines and paragraphs break automatically.
Web page addresses and e-mail addresses turn into links automatically.
Use [collapse] and [/collapse] to create collapsible text blocks. [collapse collapsed] or [collapsed] will start with the block closed.

Plain text

No HTML tags allowed.
Web page addresses and e-mail addresses turn into links automatically.
Lines and paragraphs break automatically.
Code snippets in <code>...</code> or <source>...</source> automatically will be pretty printed.
Use [collapse] and [/collapse] to create collapsible text blocks. [collapse collapsed] or [collapsed] will start with the block closed.

CAPTCHA

This question is for testing whether or not you are a human visitor and to prevent automated spam submissions.

Drag or click the correct shape to the grey "drop area".

You are here

Matlab: Plot objective and gradients (partial derivatives)

Add new comment

Filtered HTML

Plain text