Recent Posts (page 2 / 6)

by Al Danial

Solve Global Optimization Problems in MATLAB with differential_evolution()

Part 8 of the Python Is The Ultimate MATLAB Toolbox series.



Introduction

MATLAB’s fminsearch() function finds local minima of unbounded multivariable functions. With a lucky starting guess, it may even return the global minimum. For more rigirous searches of global minima, or control over the search domain, you’ll need more powerful options such as those in the Optimization Toolbox. Alternatives to the Optimization Toolbox include global optimizers from

  1. the FileExchange, specifically a differential evolution solver as used in this post, or
  2. the scipy.optimize module, provided that you’re willing to include Python in your MATLAB workflow and implement your objective function in Python.

The first option is the most convenient for MATLAB users as it provides a pure MATLAB solution. The second option may be necessary in corporate networks that restrict access to the Internet but include Python on the same computers that have MATLAB. This post covers both options.

Differential evolution

scipy.optimize has three global optimizers: shgo(), dual_annealing(), and differential_evolution(). Of the three, the most powerful is arguably differential_evolution() which implements the Storn and Price algorithm1 of the same name.

Briefly, the differential evolution algorithm works by starting with a population of random guesses to the solution. Guesses which yield better answers to the cost function are combined to form new guesses while poorly performing guesses are eliminated. The algorithm repeats this process until the guesses converge.

An optimization problem: shape alignment

Say you have two images of a triangular widget, where the widget’s position in one image is offset from the position in the other. After some image processing, perhaps with OpenCV to identify only the corners, you’re left with these two images:

Target image on the left, known image on the right.

Target image on the left, known image on the right.

Our goal is to figure out how far to translate the shape on the right until it overlaps the one on the left (we’ll ignore rotation and assume the widgets are always oriented the same way). Simple, right? Just compute the horizontal and vertical differences of a common corner such as the upper tip:

However if the widget on the left is on a noisy background and the corner identifcation code returns many spurious points such as shown here,

then problem is much harder since there’s no clear way to determine which known point corresponds to which target point.

There exist algorithms that can compute accurate offsets quickly but here we’ll take a lazy approach and let an optimizer figure this out. The cost function can be the sum of distances from each corner of our known image to the nearest point on the target image. This sum should be zero if all three corners of the known image perfectly overlay the three nearest target points.

We’ll start with three known points on a unit grid. To this we’ll add ten target points which are made of three copies of the known points plus an offset of (-0.4, -0.2) and seven additional random points. We want an optimizer to return the offset using only the locations of the original three known points and the ten target points. The code below makes a 10 x 2 array of random coordinates for the targets then arbitrarily replaces the 3rd, 6th, and 10th points with the known points plus the offset.

Python:

In : import numpy as np
In : known = np.array([[0.6, 0.5], [0.6, 0.9], [0.7,0.5]])
In : target = np.random.rand(10,2)
In : truth_offset = np.array([-0.4, -0.2])
In : target[2,:] = known[0,:] + truth_offset
In : target[5,:] = known[1,:] + truth_offset
In : target[9,:] = known[2,:] + truth_offset

MATLAB:

>> known = [0.6 0.5; 0.6 0.9; 0.7 0.5];
>> target = rand(10,2);
>> truth_offset = [-0.4 -0.2];
>> target( 3,:) = known(1,:) + truth_offset;
>> target( 6,:) = known(2,:) + truth_offset;
>> target(10,:) = known(3,:) + truth_offset;

Python and MATLAB will come up with different random target points which will make it harder to compare solutions. To keep things consistent I’ll run the MATLAB version to produce a set of values, then hard-code these into both MATLAB and Python test environments:

Python:

known = np.array([[0.6, 0.5], [0.6, 0.9], [0.7,0.5]])
target = np.array([[ 8.147236863931789e-01, 1.576130816775483e-01],
                   [ 9.057919370756192e-01, 9.705927817606157e-01],
                   [ 2.0e-01              , 3.0e-01              ],
                   [ 9.133758561390194e-01, 4.853756487228412e-01],
                   [ 6.323592462254095e-01, 8.002804688888001e-01],
                   [ 2.0e-01              , 7.0e-01              ],
                   [ 2.784982188670484e-01, 4.217612826262750e-01],
                   [ 5.468815192049838e-01, 9.157355251890671e-01],
                   [ 9.575068354342976e-01, 7.922073295595544e-01],
                   [ 3.0e-01              , 3.0e-01              ]])

MATLAB:

>> known = [0.6 0.5; 0.6 0.9; 0.7 0.5];
>> target = [ 8.147236863931789e-01, 1.576130816775483e-01;
              9.057919370756192e-01, 9.705927817606157e-01;
              2.0e-01              , 3.0e-01              ;
              9.133758561390194e-01, 4.853756487228412e-01;
              6.323592462254095e-01, 8.002804688888001e-01;
              2.0e-01              , 7.0e-01              ;
              2.784982188670484e-01, 4.217612826262750e-01;
              5.468815192049838e-01, 9.157355251890671e-01;
              9.575068354342976e-01, 7.922073295595544e-01;
              3.0e-01              , 3.0e-01              ];

Next we need to find the closest target point for each known point To do that we’ll need to compute distances from all targets to all knowns then pick the target with the smallest distance for each known point. The pairwise distance matrix computation is simple with NumPy arrays in Python and matrices in MATLAB because both support broadcasting. With broadcasting, subtracting a row vector with nR terms from a column vector with nC terms leaves us with an nR x nC matrix of all pairwise differences:

MATLAB:

>> format short e

>> known(:,1)'
   6.0000e-01   6.0000e-01   7.0000e-01

>> target(:,1)

   8.1472e-01
   9.0579e-01
   2.0000e-01
   9.1338e-01
   6.3236e-01
   2.0000e-01
   2.7850e-01
   5.4688e-01
   9.5751e-01
   3.0000e-01

then

MATLAB:

>> target(:,1) - known(:,1)'

   2.1472e-01   2.1472e-01   1.1472e-01
   3.0579e-01   3.0579e-01   2.0579e-01
  -4.0000e-01  -4.0000e-01  -5.0000e-01
   3.1338e-01   3.1338e-01   2.1338e-01
   3.2359e-02   3.2359e-02  -6.7641e-02
  -4.0000e-01  -4.0000e-01  -5.0000e-01
  -3.2150e-01  -3.2150e-01  -4.2150e-01
  -5.3118e-02  -5.3118e-02  -1.5312e-01
   3.5751e-01   3.5751e-01   2.5751e-01
  -3.0000e-01  -3.0000e-01  -4.0000e-01

The same can be done with NumPy arrays but there one must explicitly add a second dimension to one of the two column extracts with np.newaxis. Without this step, the arrays of x coordinates are merely one dimensional arrays with no concept of row or column orientation. (A detailed explanation can be found in section 11.1.12 of my book.)

Python:

In : target[:,0][:,np.newaxis] - known[:,0]
Out:
array([[ 0.21472369,  0.21472369,  0.11472369],
       [ 0.30579194,  0.30579194,  0.20579194],
       [-0.4       , -0.4       , -0.5       ],
       [ 0.31337586,  0.31337586,  0.21337586],
       [ 0.03235925,  0.03235925, -0.06764075],
       [-0.4       , -0.4       , -0.5       ],
       [-0.32150178, -0.32150178, -0.42150178],
       [-0.05311848, -0.05311848, -0.15311848],
       [ 0.35750684,  0.35750684,  0.25750684],
       [-0.3       , -0.3       , -0.4       ]])

By extending the broadcasting idea to include the delta y computation then squaring and summing, we can get the complete 10 x 3 array of squares of distances between points with just one (long) line of code:

Python:

In : dist2 = ( target[:,0][:,np.newaxis] - known[:,0] )**2 + \
   :         ( target[:,1][:,np.newaxis] - known[:,1] )**2
In : dist2
Out:
array([[0.16333506, 0.5972446 , 0.13039033],
       [0.31496628, 0.09849205, 0.26380789],
       [0.2       , 0.52      , 0.29      ],
       [0.0984183 , 0.27011778, 0.04574313],
       [0.09121548, 0.01099111, 0.09474363],
       [0.2       , 0.2       , 0.29      ],
       [0.10948469, 0.33207567, 0.18378505],
       [0.1756576 , 0.00306918, 0.1962813 ],
       [0.21319626, 0.1394304 , 0.15169489],
       [0.13      , 0.45      , 0.2       ]])

MATLAB:

>> dist2 = ( target(:,1) - known(:,1)' ).^2 + ...
           ( target(:,2) - known(:,2)' ).^2

   1.6334e-01   5.9724e-01   1.3039e-01
   3.1497e-01   9.8492e-02   2.6381e-01
   2.0000e-01   5.2000e-01   2.9000e-01
   9.8418e-02   2.7012e-01   4.5743e-02
   9.1215e-02   1.0991e-02   9.4744e-02
   2.0000e-01   2.0000e-01   2.9000e-01
   1.0948e-01   3.3208e-01   1.8379e-01
   1.7566e-01   3.0692e-03   1.9628e-01
   2.1320e-01   1.3943e-01   1.5169e-01
   1.3000e-01   4.5000e-01   2.0000e-01

meaning, for example, the square of the distance between the 8th target point (== 8th row) and 3rd known point (== 3rd column) is 1.9628e-01. All that remains is to take the square root of every term.

Setting up the cost function

Now that we know how far apart each target point is from each known point, we just need to find the closest ones then sum their distances. Taking the minumum over each column, then summing gives us the cost function:

Python:

In : dist  = np.sqrt( ( target[:,0][:,np.newaxis] - known[:,0] )**2 +
   :                  ( target[:,1][:,np.newaxis] - known[:,1] )**2 )

In : np.min(dist, axis=0)
Out: array([0.30201901, 0.05540018, 0.21387643])

In : np.min(dist, axis=0).sum()
Out: 0.5712956164197415

MATLAB:

>> dist = sqrt( ( target(:,1) - known(:,1)' ).^2 + ...
                ( target(:,2) - known(:,2)' ).^2 );

>> min(dist)

    0.3020    0.0554    0.2139

>> sum(min(dist))

    0.5713

Optimizers need to call the function that’s to be minimized so we’ll rewrite the expressions above as functions called cost()2. We merely need to wrap the distance summation computation into a function:

Python:

import numpy as np
def cost(known, target):
    dist = np.sqrt( ( target[:,0][:,np.newaxis] - known[:,0] )**2 + 
                    ( target[:,1][:,np.newaxis] - known[:,1] )**2 )
    return np.min(dist, axis=0).sum()

MATLAB:

function [sum_dist] = cost(known, target)
    dist = sqrt( ( target(:,1) - known(:,1)' ).^2 + ...
                 ( target(:,2) - known(:,2)' ).^2 );
    sum_dist = sum(min(dist));
end

Visualizing the cost function

It’s instructive to see how the cost function varies across the search domain but this is only practical if the cost function has two or fewer independent variables and is inexpensive to compute. The cost function chosen for this post was obviously tailored to meet these criteria. I’d have liked to allow the target shape to be oriented at an arbitrary angle but then the offset would consist of three independent variables, x, y, and phi, and the cost function would need to be rendered as a 3D cube. That’s possible but requires an advanced visualization tool that supports cutting planes.

Here are functions to display a contour plot of the cost function across the search domain. Both functions descretize the search domain of -1 <= x <= 1 and -1 <= y <= 1 into an evenly spaced grid of 100 x 100 points then compute the cost function at each of these. Afterwards we’ll view the 100 x 100 cost terms on a contour plot.

Python:

import numpy as np
import matplotlib.pyplot as plt
def plot_cost(known, target): # {{{
    N = 100
    [offset_x, offset_y] = np.meshgrid(np.linspace(-1,1,N),np.linspace(-1,1,N))
    C = np.zeros_like(offset_x) # cost at each offset point
    offset_target = np.zeros_like(target)
    for r in range(N):
        for c in range(N):
            offset = np.array([offset_x[r,c], offset_y[r,c]])
            offset_target = target + offset
            C[r,c] = cost(known, offset_target)
    fig, ax = plt.subplots()
    plt.contourf(offset_x, offset_y, C, levels=20)
    ax.set_aspect('equal')
    ax.set_xlim([-1,1.0])
    ax.set_xlabel('X offset')
    ax.set_ylabel('Y offset')
    ax.set_ylim([-1,1.0])
    plt.colorbar()
    plt.show()

MATLAB:

function plot_cost(known, target)
    N = 100;
    [offset_x, offset_y] = meshgrid(linspace(-1,1,N), linspace(-1,1,N));
    C = zeros(size(offset_x)); % cost at each offset point
    offset_target = zeros(size(target));
    for r = 1:N
        for c = 1:N
            offset = [offset_x(r,c), offset_y(r,c)];
            offset_target = target + offset;
            C(r,c) = cost(known, offset_target);
        end
    end
    levels = 20;
    contourf(offset_x, offset_y, C, levels)
    colorbar()
end

For our first plot let’s start with the easy case with just three target points which are offset from three known points:

Python:

import numpy as np
import matplotlib.pyplot as plt
known = np.array([[0.6, 0.5], [0.6, 0.9], [0.7,0.5]])
truth_offset = np.array([-0.4, -0.2])
target = known + truth_offset;
plot_cost(known, target)

Cost function for three targets, Python.

Cost function for three targets, Python.

There’s a fairly distinct minimum around (0.4, 0.2). This is the amount we’d need to add to our target points to make them overlay the known points. Looking at it from the other direction, (-0.4, -0.2) was added to our known points to produce the targets. The cost function is fairly smooth with a clear minimum and should not present much challenge for gradient-based optimizers such as MATLAB’s fminsearch().

Now let’s plot the harder case of 10 target points—3 from shifting the known points, plus the 7 noise points:

Python:

import numpy as np
import matplotlib.pyplot as plt
known = np.array([[0.6, 0.5], [0.6, 0.9], [0.7,0.5]])
target = np.array([[ 8.147236863931789e-01, 1.576130816775483e-01],
                   [ 9.057919370756192e-01, 9.705927817606157e-01],
                   [ 2.0e-01              , 3.0e-01              ],
                   [ 9.133758561390194e-01, 4.853756487228412e-01],
                   [ 6.323592462254095e-01, 8.002804688888001e-01],
                   [ 2.0e-01              , 7.0e-01              ],
                   [ 2.784982188670484e-01, 4.217612826262750e-01],
                   [ 5.468815192049838e-01, 9.157355251890671e-01],
                   [ 9.575068354342976e-01, 7.922073295595544e-01],
                   [ 3.0e-01              , 3.0e-01              ]])
plot_cost(known, target)

Cost function for ten targets, Python.

Cost function for ten targets, Python.

This time the global minimum is hard to identify. If we add even more noise points the plot will have more local minima and the optimization problem becomes even more challenging.

For completeness here are the MATLAB versions of the two plots:

MATLAB:

known = [0.6 0.5; 0.6 0.9; 0.7 0.5];
truth_offset = [-.4 -.2];
target3 = known + truth_offset;
plot_cost(known, target3)

Cost function for three targets, MATLAB.

Cost function for three targets, MATLAB.

MATLAB:

known = [0.6 0.5; 0.6 0.9; 0.7 0.5];
target = [ 8.147236863931789e-01, 1.576130816775483e-01;
           9.057919370756192e-01, 9.705927817606157e-01;
           2.0e-01              , 3.0e-01              ;
           9.133758561390194e-01, 4.853756487228412e-01;
           6.323592462254095e-01, 8.002804688888001e-01;
           2.0e-01              , 7.0e-01              ;
           2.784982188670484e-01, 4.217612826262750e-01;
           5.468815192049838e-01, 9.157355251890671e-01;
           9.575068354342976e-01, 7.922073295595544e-01;
           3.0e-01              , 3.0e-01              ];
plot_cost(known, target)

Cost function for ten targets, MATLAB.

Cost function for ten targets, MATLAB.

Finding the global minimum

We’re ready to put our optimizers to work. Before we jump to differential evolution, let’s try the fminsearch() function in MATLAB, first with the three target points:

MATLAB:

>> start_offset = [0 0];
>> fn = @(offset) cost(known, target3 + offset);
>> fminsearch(fn, start_offset)

    0.4000    0.2000

Looks good. Now the much harder case with the ten target points:

MATLAB:

>> start_offset = [0 0];
>> fn = @(offset) cost(known, target + offset);
>> fminsearch(fn, start_offset)

    0.0285   -0.2725

Not so good—(0.0285, -0.2725) is one of the several local minima. The Optimization Toolbox has functions to find the global minimum but if that’s not available, a good alternative is a differential evolution based solver from the FileExchange or from scipy.optimize.

Python solution with scipy.optimize.differential_evolution()

SciPy’s differential_evolution() calling arguments are similar to MATLAB’s fminsearch() but include an additional array of bounds on the independent variables (the x and y values of the offset). The first argument is a handle to a function, we’ll call it fn() to match the @fn() anonymous function used in MATLAB since it serves the same purpose—in both Python and MATLAB, fn calls our cost function with trial values for independent variables plus other required arguments, namely the known and target arrays, passed through with the keyword argument args. Here’s what that looks like:

Python:

def fn(independent_vars, *data):
    x_offset, y_offset = independent_vars
    known, target = data
    offset_target = target + np.array([x_offset, y_offset])
    return cost(known, offset_target)

bounds=([-1.0, 1.0], [ -1.0, 1.0])

result = differential_evolution(fn, bounds, args=(known, target))

Solutions for the three and ten point target cases return the correct global minimum for both cases:

Python:

In : differential_evolution(fn, bounds, args=(known, target3))
Out:
     fun: 4.052626611931048e-16
 message: 'Optimization terminated successfully.'
    nfev: 3033
     nit: 98
 success: True
       x: array([0.4, 0.2])

In : differential_evolution(fn, bounds, args=(known, target))
Out:
     fun: 3.3664710476199897e-16
 message: 'Optimization terminated successfully.'
    nfev: 2913
     nit: 94
 success: True
       x: array([0.4, 0.2])

Differential evolution involves random guesses for the independent variables so the number of function evaluations (nfev in the solution structure) varies between runs. Oddly, in this instance, the easier case of three target points required more evaluations than the harder case of ten points. If the cost function is expensive—for example a structural analysis run using Ansys or NASTRAN to find the maximum stress in a finite element model—you may want to modify the differential_evolution() optional parameters to find a suitable solution with fewer evaluations, or to run on multiple processors.

Individual elements of the returned structure are accessed through dot notation, for example,

Python:

In : result = differential_evolution(fn, bounds, args=(known, target3))
In : result.nfev
Out: 3213

In : result.x
Out: array([0.4, 0.2])

The complete Python program with all code above is optim_w_DE.py.

MATLAB solution with differential_evolution() from the FileExchange

Markus Buehren posted a thorough implementation of differential evolution entirely in MATLAB on the FileExchange. I didn’t spend enough time with it to learn how to pass arbitrary inputs to the cost function so the implementation below embeds the known and target arrays directly in cost().

I expanded the zip file to the directory /home/al/matlab/DE on my computer. To use Markus’s solver this directory must be in the MATLAB search path so the script below begins with an addpath() to this location.

MATLAB: optim_w_fileexchange.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
% file: optim_w_fileexchange.m
addpath('/home/al/matlab/DE'); % location of unzipped files
optimInfo.title = 'Shape alignment';
tolerance = 0.0001;
paramDefCell = {
	'parameter1', [-1 1], tolerance  % offset_x
	'parameter2', [-1 1], tolerance  % offset_y
};
objFctParams.offset_x = 0;
objFctParams.offset_y = 0;
objFctHandle = @cost;
DEParams = getdefaultparams;
DEParams.maxiter = 50;
objFctSettings = 100;
emailParams = [];  % no email
[bestmem, bestval, bestFctParams, nrOfIterations, resultFileName] = ...
    differentialevolution(DEParams, paramDefCell, objFctHandle,     ...
    objFctSettings, objFctParams, emailParams, optimInfo);

function [sum_dist] = cost(scale, params)
    offset_x = params.parameter1(1);
    offset_y = params.parameter2(1);
    known = [0.6 0.5; 0.6 0.9; 0.7 0.5];
    target = [ 8.147236863931789e-01, 1.576130816775483e-01;
               9.057919370756192e-01, 9.705927817606157e-01;
               2.0e-01              , 3.0e-01              ;
               9.133758561390194e-01, 4.853756487228412e-01;
               6.323592462254095e-01, 8.002804688888001e-01;
               2.0e-01              , 7.0e-01              ;
               2.784982188670484e-01, 4.217612826262750e-01;
               5.468815192049838e-01, 9.157355251890671e-01;
               9.575068354342976e-01, 7.922073295595544e-01;
               3.0e-01              , 3.0e-01              ];
    offset_target = target + [offset_x, offset_y];
    dist = sqrt( ( offset_target(:,1) - known(:,1)' ).^2 + ...
                 ( offset_target(:,2) - known(:,2)' ).^2 );
    sum_dist = scale*sum(min(dist));
end

Running this script yields the correct solution

MATLAB:

                                   :
'Shape alignment' finished after given maximum number of 50 iterations.
 
Elapsed time:                   0.37 seconds
Number of function evaluations: 961
Mean function evaluation time:  0.0003928 seconds

Final result after 50 iteration(s).
Best member:
parameter1 = 0.4; 
parameter2 = 0.2; 
Evaluation value: 5.10014e-14

Results are saved in file cost_result_01.mat.

and additionally brings up plots and writes .mat files to the working directory.

MATLAB solution with scipy.optimize.differential_evolution()

If the FileExchange implementation of differential evolution isn’t an option for you, use the implementation in SciPy instead. Although scipy.optimize.differential_evolution() can’t be called from MATLAB directly—it requires a Python evaluation function, the fn() mentioned above—the call can be put in a Python wrapper and the wrapper called from MATLAB. This means, however, the cost function must also be implemented in Python (or I suppose could use the MATLAB engine in Python to call the MATLAB cost function from Python although I’ve never attempted this). Here’s what the wrapper looks like for our problem:

Python: diffev_bridge.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# file: diffev_bridge.py
import numpy as np
from scipy.optimize import differential_evolution

def cost(known, target):
    dist = np.sqrt( ( target[:,0][:,np.newaxis] - known[:,0] )**2 +
                    ( target[:,1][:,np.newaxis] - known[:,1] )**2 )
    return np.min(dist, axis=0).sum()

def fn(independent_vars, *data):
    x_offset, y_offset = independent_vars
    known, target = data
    offset_target = target + np.array([x_offset, y_offset])
    return cost(known, offset_target)

def call_DE(bounds, known, target):
    args = (known, target)
    DE_result = differential_evolution(fn, bounds, args=args)
    soln = { 'success' : DE_result.success,
             'message' : DE_result.message,
             'x' : DE_result.x,
             'nfev' : DE_result.nfev,
           }
    return soln

The MATLAB invocation to this wrapper is

MATLAB (>= 2020b): optim_w_python.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
np = py.importlib.import_module('numpy');
diffev = py.importlib.import_module('diffev_bridge');

known = [0.6 0.5; 0.6 0.9; 0.7 0.5];
offset = [-.4 -.2];
target = [ 8.147236863931789e-01, 1.576130816775483e-01;
           9.057919370756192e-01, 9.705927817606157e-01;
           2.0e-01              , 3.0e-01              ;
           9.133758561390194e-01, 4.853756487228412e-01;
           6.323592462254095e-01, 8.002804688888001e-01;
           2.0e-01              , 7.0e-01              ;
           2.784982188670484e-01, 4.217612826262750e-01;
           5.468815192049838e-01, 9.157355251890671e-01;
           9.575068354342976e-01, 7.922073295595544e-01;
           3.0e-01              , 3.0e-01              ];
target3 = known + offset;
bounds = py.tuple({[-1, 1], [-1, 1]});
result = diffev.call_DE(bounds, np.array(known), np.array(target3))
result = diffev.call_DE(bounds, np.array(known), np.array(target))

running it returns a solution dictionary for the three point and ten point cases:

MATLAB (>= 2020b):

>> optim_w_python
result = 
  Python dict with no properties.
    {'success': True, 
     'message': 'Optimization terminated successfully.', 
     'x': array([0.4, 0.2]), 
     'nfev': 3153}

result = 
  Python dict with no properties.
    {'success': True, 
     'message': 'Optimization terminated successfully.', 
     'x': array([0.4, 0.2]), 
     'nfev': 3033}

Individual items from the solution dictionary can be accessed using the dictionary .get() method. Alternatively, convert the dictionary to a native MATLAB struct with py2mat.m:

MATLAB (>= 2020b):

>> offset = double(result.get('x'))

    0.4000    0.2000

>> soln = py2mat(result)

  struct with fields:

    success: 1
    message: "Optimization terminated successfully."
          x: [0.4000 0.2000]
       nfev: 2943

>> soln.x

    0.4000    0.2000

Avoiding MATLAB crashes on Linux

MATLAB 2020b on Ubuntu 20.04 may crash at diffev.call_DE() when running optim_w_python.m. The solution is to start MATLAB from a suitably-configured Python environment and to include these lines in startup.m:

MATLAB (>= 2020b):

py.sys.setdlopenflags(int32(10));
os = py.importlib.import_module('os');
os.environ.update(pyargs('KMP_DUPLICATE_LIB_OK','TRUE'))

Additional Python optimization packages

Optimization is a massive topic that spans problem classes that differential evolution cannot solve. See Keivan Tafakkori’s list of optimization packages for Python for an overview of additional solvers.

Join me again on September 17, 2022 to see how to solve simultaneous nonlinear equations in MATLAB using SciPy’s scipy.optimize.root() function.


  1. Storn, R.; Price, K. (1997). “Differential evolution - a simple and efficient heuristic for global optimization over continuous spaces”. Journal of Global Optimization. 11 (4): 341–359. doi:10.1023/A:1008202821328. S2CID 5297867. ↩︎

  2. If you are using the Audio Toolbox, Communications Toolbox, or DSP System Toolbox you’ll want to call the function something other than cost() since these toolboxes also have a cost() function. ↩︎

by Al Danial

Work with Redis in MATLAB Using the redis-py Python Module

Part 7 of the Python Is The Ultimate MATLAB Toolbox series.

Redis, the “remote dictionary server”, is a fast network-based in-memory database for key/value pairs. In addition to high performance and ability to store structured values, Redis implements a publish/subscribe service so that applications can be notified of changes to keys they subscribe to. The MATLAB Production Server has a Redis interface. Python does too—several, actually. In this article I’ll demonstrate a MATLAB interface to Redis via the the redis-py Python module.

Why would a MATLAB user care about Redis?

Redis may prove useful to you if you have a MATLAB application that needs to share values rapidly with other programs on your network. The programs involved use Redis as a remote memory bank for shared variables. Interaction with a Redis server happens with the Redis communication API which is available for more than 50 languages including Python—and therefore, by extension, MATLAB as well. While that sounds complicated to set up, it is often easier than writing other networking code, whether directly with TCP/IP sockets, or with a networking framework such as MPI. It is certainly cleaner and faster than communicating through the file system.

Note: Redis stores all values as character (or byte) arrays. For this reason it is poorly-suited for exchanging numeric arrays as might be needed by a computationally intensive parallel MATLAB program. Stick with MPI from the Parallel Computing Toolbox for such applications (or wait for my Nov. 12, 2022 article on using Dask with MATLAB).

If your MATLAB projects work in isolation, that is, they don’t need to get or provide information to other applications, chances are you won’t need a network-based key/value store. If that’s the case, the rest of this article may not interest you.

Install redis-py

You’ll need the redis-py Python module in your installation to call it from MATLAB.

NOTE: There’s a confusing overlap of names for the Python Redis package, the Redis application, and Anaconda’s Redis bundle.

Is it redis or redis-py ?

  • The primary Python Redis module is developed at https://github.com/redis/redis-py . Although the name on Github is redis-py, once you install this module, you import it with import redis (not import redis-py)

  • In Anaconda distributions, this Python module is installed with conda install redis-py (not conda install redis)

  • In Anaconda distributions, the full-up Redis application including the server and command line tools can be installed with conda install redis

  • In other Python distributions, the module can be installed with python -m pip install redis (not python -m pip install redis-py)

Check your installation

Check your Python installation to see if it has the redis-py Python module by starting the Python command line REPL and importing it—with import redis, not import redis-py. If it imports successfully you can also check the version number.

> python
>>> import redis
>>> redis.__version__
'3.5.3'

If you get a ModuleNotFoundError, you’ll need to install it using one of the methods listed above.

Start a Redis server

Skip this section if you already have a Redis server running somewhere on your network (or on your own computer).

To begin, install the Redis application by following the Redis installation instructions, using your operating system’s package manager (for example on Ubuntu you’d need apt install redis-tools redis-server) or, if you have the Anaconda Python installation, by installing Redis with conda:

> conda install redis

This provides both the server, redis-server and command line tools such as the redis-cli client program. Next, open a terminal and start the server:

> redis-server  redis.conf --port 12987

The configuration file redis.conf and port value --port 12987 are optional. If you don’t specify them the server will use a default configuration and run on port 6379. You’ll need to create a config file if you want to add password protection and persist the memory database to disk, among many other options.

If all is well you’ll see some ASCII art:

                _._
           _.-``__ ''-._
      _.-``    `.  `_.  ''-._           Redis 5.0.3 (00000000/0) 64 bit
  .-`` .-```.  ```\/    _.,_ ''-._                                   
 (    '      ,       .-`  | `,    )     Running in standalone mode
 |`-._`-...-` __...-.``-._|'` _.-'|     Port: 12987
 |    `-._   `._    /     _.-'    |     PID: 2453308
  `-._    `-._  `-./  _.-'    _.-'                                   
 |`-._`-._    `-.__.-'    _.-'_.-'|
 |    `-._`-._        _.-'_.-'    |           http://redis.io 
  `-._    `-._`-.__.-'_.-'    _.-'                                   
 |`-._`-._    `-.__.-'    _.-'_.-'|
 |    `-._`-._        _.-'_.-'    |
  `-._    `-._`-.__.-'_.-'    _.-'                                   
      `-._    `-.__.-'    _.-'                                       
          `-._        _.-'                                           
              `-.__.-'                                               

2453308:M 30 Jul 2022 14:30:20.334 # Server initialized
2453308:M 30 Jul 2022 14:30:20.334 * Ready to accept connections

Connect to a Redis server

We’ll use the redis-cli command line application to test our ability to interact with the Redis server. It starts a read-evaluate-print loop (REPL) that accepts Redis commands:

(matpy) » redis-cli -h localhost -p 12987
localhost:12987>

The prompt above indicates a successful connection to the server running on the same computer. Had the connection failed, perhaps because the server isn’t running, the prompt would look like this:

(matpy) » redis-cli -h localhost -p 12987
Could not connect to Redis at localhost:12987: Connection refused
not connected>

The -h localhost switch is redundant as the default host is the current computer. It appears here explicitly to give a sense of how one would connect to a Redis server on a different computer. Also the -p 12987 can be omitted if your Redis server is using the default port (6379).

Add or update a key/value pair

Within a successfully-connected redis-cli session we can add new keys and values with the SET command:

localhost:12987> set sensor_T712_temperature 23.4
OK
localhost:12987>

Any Redis-capable application on our network can now retrieve the value for sensor_T712_temperature. To update the value, simply set it again:

localhost:12987> set sensor_T712_temperature 37.9
OK
localhost:12987>

Get a key’s value

We can retrieve a key’s value with the GET command

localhost:12987> get sensor_T712_temperature
"37.9"
localhost:12987>

Note that the return value is a character array, not a floating point number.

Batch redis-cli commands

Exit the redis-cli REPL session above with <ctrl>-c, <ctrl>-d, or by enterring either exit or quit. The command line client program accepts complete Redis commands without needing to entering the REPL. This makes it possible to write Redis-capable shell or Windows .bat scripts, or in fact be used by any programming language capable of making a system call. We’ll start with the ping command which merely tests the connection. A reply of PONG indicates success:

> redis-cli -h localhost -p 12987 ping
PONG

> redis-cli -h localhost -p 12987 set sensor_T712_temperature 40.2
OK

> redis-cli -h localhost -p 12987 get sensor_T712_temperature
"40.2"

Now that we’ve checked that our Redis server is up and we can interact with it, let’s do the interactions programmatically:

Python: set_get.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
#!/usr/bin/env python
# set_get.py
import redis
R = redis.Redis(host='localhost',port=12987,decode_responses=True)

try:
    R.ping()
except redis.exceptions.ConnectionError as e:
    print(f'Is Redis running?  Unable to connect {e}')
    sys.exit(1)
print(f'Connected to server.')

R.set('sensor_T712_temperature', 39.8)
retrieved_temp = R.get('sensor_T712_temperature')
print(f'sensor_T712_temperature = {retrieved_temp}')
print(f'(as a number)           = {float(retrieved_temp):.3f}')

A note about the decode_responses=True option: without it, programmatic reads of Redis keys and values receive byte arrays. In most cases you’d rather have strings, not byte arrays. You can get a string from a byte array by applying the .decode() method to the byte array, but including the decode_responses=True in the initial connection call spares you from having to add .decode() to every Redis return variable.

In any event, you still have to explicitly typecast the string temperature value if you want to use it numerically.

Running the Python program gives:

> ./set_get.py
Connected to server.
sensor_T712_temperature = 39.8
(as a number)           = 39.800

Now in MATLAB:

MATLAB: set_get.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
% set_get.m
redis = py.importlib.import_module('redis');
R = redis.Redis(pyargs('host','localhost','port',int64(12987),...
                       'decode_responses',py.True));
try
  R.ping();
  fprintf('Connected to server.\n')
catch EO
  fprintf('Is Redis running?  Unable to connect: %s\n', EO.message)
end

R.set('sensor_T712_temperature', 39.8);
retrieved_temp = R.get('sensor_T712_temperature');
fprintf('sensor_T712_temperature = %s\n', string(retrieved_temp))
as_number = double(string(retrieved_temp));
fprintf('(as a number)           = %6.3f\n', as_number)

>> set_get
Connected to server.
sensor_T712_temperature = 39.8
(as a number)           = 39.800

Set/get structured data

The ‘value’ part of a Redis key/value pair may be a structure. Supported structures are a string, set, sorted set, list, hash, bitmap, bitfield, hyperloglog, stream, and geospatial index. The example below demonstrates setting and getting a hash since conceptually this structure resembles a MATLAB struct.

Add a hash from the command line

The Redis command HSET adds a hash. This command

> redis-cli -h localhost -p 12987 hset newmark alpha 0.25 beta 0.5 gamma 'is unset'

adds a key newmark with fields alpha, beta, and gamma corresponding to a struct like this in MATLAB:

>> newmark
  struct with fields:
    alpha: 0.2500
     beta: 0.5000
    gamma: "is unset"

Bear in mind though the numeric values in the MATLAB struct are saved as strings in the Redis hash. The Redis hash can be retrieved in its entirety with HGETALL

> redis-cli -h localhost -p 12987 hgetall newmark
1) "alpha"
2) "0.25"
3) "beta"
4) "0.5"
5) "gamma"
6) "is unset"

or selectively with HMGET (where M is for “multi-values”) by providing a list of the desired fields:

> redis-cli -h localhost -p 12987 hmget newmark alpha gamma
1) "0.25"
2) "is unset"

Programmatically, the .hgetall() function returns a Python dictionary in both Python and MATLAB. This is convenient when the values are strings but not that helpful if the values have to be typecast to numbers. In that case it’s an equal amount of work to instead calling the .hmget() method to just get back a list of values (once again, byte arrays) and convert those.

Python: get_hash.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
#!/usr/bin/env python
# get_hash.py
import redis
R = redis.Redis(host='localhost',port=12987,decode_responses=True)

try:
    R.ping()
except redis.exceptions.ConnectionError as e:
    print(f'Is Redis running?  Unable to connect {e}')
    sys.exit(1)
print(f'Connected to server.')

values = R.hmget('newmark', ['alpha', 'beta', 'gamma'])
newmark = { 'alpha' : float(values[0]),
            'beta'  : float(values[1]),
            'gamma' : values[2] }
print(newmark)

> get_hash.py
Connected to server.
{'alpha': 0.25, 'beta': 0.5, 'gamma': 'is unset'}

MATLAB: get_hash.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
% get_hash.m
redis = py.importlib.import_module('redis');
R = redis.Redis(pyargs('host','localhost','port',int64(12987),...
                       'decode_responses',py.True));
try
  R.ping();
  fprintf('Connected to server.\n')
catch EO
  fprintf('Is Redis running?  Unable to connect: %s\n', EO.message)
end

values = R.hmget('newmark', {'alpha', 'beta', 'gamma'});
newmark = struct;
newmark.alpha = double(string(values{1}));
newmark.beta  = double(string(values{2}));
newmark.gamma = string(values{3});
disp(newmark)

>> get_hash
Connected to server.
    alpha: 0.2500
     beta: 0.5000
    gamma: "is unset"

Set/get performance

The following programs set, then get 100,000 keys with integer values from 0 to 99,999. The time to .decode() the returned byte array is also included for the get time since this operation is nearly always necessary. The MATLAB code is more than 2x slower on sets and a baffling 5x slower on gets.

Python: set_get_speed.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
#!/usr/bin/env python
# set_get_speed.py
import redis
import time
R = redis.Redis(host='localhost',port=12987)

try:
    R.ping()
except redis.exceptions.ConnectionError as e:
    print(f'Is Redis running?  Unable to connect {e}')
    sys.exit(1)
print(f'Connected to server.')

N = 100_000

T_s = time.time()
for i in range(N):
    R.set(f'K_%06{i}', i)
T_e = time.time()
print(f'{N} sets for integer values took {T_e-T_s:.3f} s')

T_s = time.time()
for i in range(N):
    v = R.get(f'K_%06{i}')
    v = v.decode()
T_e = time.time()
print(f'{N} gets for integer values took {T_e-T_s:.3f} s')

MATLAB: set_get_speed.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
% set_get_speed.m
redis = py.importlib.import_module('redis');
R = redis.Redis(pyargs('host','localhost','port',int64(12987),...
                       'decode_responses',py.True));
try
  R.ping();
  fprintf('Connected to server.\n')
catch EO
  fprintf('Is Redis running?  Unable to connect: %s\n', EO.message)
end

N = 100000;

tic
for i = 1:N
    key = sprintf('K_%06d', i-1);
    R.set(key, i-1);
end
fprintf('%d sets for integer values took %.3f s\n', N, toc)

tic
for i = 1:N
    key = sprintf('K_%06d', i-1);
    v = R.get(key);
end
fprintf('%d gets for integer values took %.3f s\n', N, toc)

Python:

> ./set_get_speed.py
Connected to server.
100000 sets for integer values took 8.377 s
100000 gets for integer values took 7.973 s

MATLAB:

>> set_get_speed
Connected to server.
100000 sets for integer values took 21.580 s
100000 gets for integer values took 44.525 s

Despite the lower performance, a MATLAB Redis set followed by a get takes less than a millisecond.

Subscribe to key changes

A Redis power-feature is its support for subscriptions. When a program subscribes to a key change, Redis notifies the program each time the key’s value changes. The first step involves notifying Redis that you wish to monitor “keyspace events”, namely, changes to keys as well as key expiration and removal. The most common type of keyspace configuration option is “KEA”. It is set from the command line like so:

redis-cli -h localhost -p 12987 config set notify-keyspace-events KEA

or programmatically like this in Python:

import redis
R = redis.Redis(host='localhost',port=12987,decode_responses=True)
R.config_set('notify-keyspace-events', 'KEA')

or this MATLAB:

redis = py.importlib.import_module('redis');
R = redis.Redis(pyargs('host','localhost','port',int64(12987),...
                       'decode_responses',py.True));
R.config_set('notify-keyspace-events', 'KEA')

The next step is to define a string pattern that matches one or more keys to watch for. For this example, say our Redis instance has a collection of keys storing temperature, humidity, and light data from an array of sensors. Their names might be

sensor_T712_temperature
sensor_T714_temperature
sensor_T716_illumination
sensor_T712_humidity

To have our application to act on changes in any of the temperature keys, our subscription would look like this:

Sub = R.pubsub()
Sub.psubscribe('__keyspace@0__:sensor_*_temperature')

After setting up the subscription, the code waits for notifications. The most direct way of doing this is in a loop that sleeps (or does something else) when no notifications come in:

Python: subscribe.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#!/usr/bin/env python
# subscribe.py
import time
import sys
import redis
R = redis.Redis(host='localhost',port=12987,decode_responses=True)

try:
    R.ping()
except redis.exceptions.ConnectionError as e:
    print(f'Is Redis running?  Unable to connect {e}')
    sys.exit(1)
print(f'Connected to server.')

# modify configuration to enable keyspace notification
R.config_set('notify-keyspace-events', 'KEA')

Sub = R.pubsub()
keys_to_match = 'sensor_*_temperature'
Sub.psubscribe('__keyspace@0__:' + keys_to_match)
while True:
    try:
        message = Sub.get_message()
    except redis.exceptions.ConnectionError:
        print('lost connection to server')
        sys.exit(1)
    if message is None:
        time.sleep(0.1)
        continue
    keyname = message['channel'].replace('__keyspace@0__:','')
    if keyname == keys_to_match:
        # initial subscription value
        continue
    # gets here if any of the matching keys was updated
    value = R.get(keyname)
    print(f'{keyname} = {value}')

MATLAB: subscribe.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
% subscribe.m
redis = py.importlib.import_module('redis');
R = redis.Redis(pyargs('host','localhost','port',int64(12987),...
                       'decode_responses',py.True));
try
  R.ping();
  fprintf('Connected to server.\n')
catch EO
  fprintf('Is Redis running?  Unable to connect: %s\n', EO.message)
end

% modify configuration to enable keyspace notification
R.config_set('notify-keyspace-events', 'KEA');

Sub = R.pubsub();
keys_to_match = "sensor_*_temperature";
Sub.psubscribe("__keyspace@0__:" + keys_to_match);

while 1
   try
       message = Sub.get_message();
   catch EO
       fprintf('lost connection to server: %s\n', EO.message)
       break
   end
   if message == py.None
       pause(0.1)
       continue
   end
   keyname = replace(message{'channel'}, '__keyspace@0__:','');
   if keyname == keys_to_match
       % initial subscription value
       continue
   end
   % gets here if any of the matching keys was updated
   value = string(R.get(keyname));
   fprintf('%s = %s\n', string(keyname), value)
end

When these programs run, they spend most of their time sleeping (line 27 in both programs for time.sleep(0.1) and pause(0.1)). Of course the sleeps should be replaced by calls to actual work if your application has better things to do. Then, when the clients receive a notification from Redis that any of the temperature sensors was updated (line 22 in Python and line 21 in MATLAB), the message variable is populated with the new temperature and the print statements at lines 35 (Python) and 37 (MATLAB) are executed.

The actual duration of the sleeps—or the other work done while not listening for the key changes—is not important. Keyspace change notifications are buffered and the clients will still be informed even if the clients call message = Sub.get_message() long after the keys actually changed.

We can run the Python and MATLAB programs simultaneously and watch them react to command line updates of some keys. A MATLAB session responsing to command line key sets such as

> redis-cli -p 12987
127.0.0.1:12987> set sensor_A_temperature 101.5
OK
127.0.0.1:12987> set sensor_B_temperature 101.6
OK
127.0.0.1:12987> set sensor_B_humidity 44.3
OK
127.0.0.1:12987> set sensor_C_temperature 101.7
OK

looks like this

>> subscribe
Connected to server.
sensor_A_temperature = 101.5
sensor_B_temperature = 101.6
sensor_C_temperature = 101.7

Note the program did not respond to a change to sensor_B_humidity since this key name doesn’t match our pattern of sensor_*_temperature.

Join me again on September 3, 2022 to see how MATLAB can take advantage of SciPy’s powerful “differential evolution” optimizer, scipy.optimize.differential_evolution.