Notebook

Brute Force Optimization (brute2fine)¶

One of the main applications of dmipy is the ability to fit any model combination to measured dMRI data. To find model parameters that appropriately describe the underlying tissue configuration, it is important to have optimization algorithms that are able to find the globally best solution based on the given data.

The most straight-forward implementation of global optimization in dmipy is brute-force. There are two implementations:

GlobalBruteOptimizer precomputes a global parameter and accompanying signal grid that is re-used when fitting different voxels. It only returns the nearest precomputed parameter combination based on the sum-squared difference between simulated signal grid and the measured data. For fitting many voxels without any prior initial parameter guess this is the fastest way to get a good intial guess. The result can then be passed to the following Brute2FineOptimizer that can refine the solution to a local minimum using gradient-based methods.
Brute2FineOptimizer does a separate brute-force optimization for every voxel and then refines the result to a local minimum using gradient-based methods. This approach is slower, but allows for (possibly partial) voxel-varying initial parameter guesses to be passed to the optimizer. Brute2FineOptimizer then only does brute-force optimization on the parameters that did not get an initial guess, and then refines the solution.

To illustrate how these optimizers work we will first set up a simple MultiCompartmentMicrostructureModel and load the hcp acquisition scheme as example

In [1]:

from dmipy.core.modeling_framework import MultiCompartmentModel
from dmipy.signal_models import gaussian_models
from dmipy.data.saved_acquisition_schemes import wu_minn_hcp_acquisition_scheme
zeppelin = gaussian_models.G2Zeppelin()
mc_model = MultiCompartmentModel([zeppelin])
scheme = wu_minn_hcp_acquisition_scheme()

GlobalBruteOptimizer using precomputed signal grid¶

Normally, these optimizers are only called inside the mc_model.fit() function, but we can also call them separately. The GlobalBruteOptimizer has a very straightforward design: given a model and an acquisition scheme, it will equally sample every model parameter between the minimum and maximum parameter ranges in 'Ns' steps. However, it treats orientation parameter 'mu' differently, where it uses a spherical grid of 'N_sphere_samples' points.

In [2]:

from dmipy.optimizers import brute2fine
Ns = 5  # number of equal sampling points for standard parameters
N_sphere_samples = 30  # number of orientations that are sampled for orientation 'mu'
global_brute_optimizer = brute2fine.GlobalBruteOptimizer(mc_model, scheme,
                                                         Ns=Ns, N_sphere_samples=N_sphere_samples)

We can see it computed a parameter and signal grid for 750 different parameter combinations

In [3]:

print 'parameter grid size:', global_brute_optimizer.parameter_grid.shape
print 'signal grid size:', global_brute_optimizer.signal_grid.shape

parameter grid size: (750, 4)
signal grid size: (750, 288)

Showing mc_model's parameters, we can see that this corresponds to $5 * 5 * 30 = 750$ samples for lambda_par, lambda_perp and mu.

In [4]:

mc_model.parameter_cardinality

Out[4]:

OrderedDict([('G2Zeppelin_1_lambda_perp', 1),
             ('G2Zeppelin_1_mu', 2),
             ('G2Zeppelin_1_lambda_par', 1)])

If we just take one of the simulated signal grid points as input for the optimizer, we can see that it will return the corresponding parameter combination.

In [5]:

test_signal = global_brute_optimizer.signal_grid[710]
test_parameters = global_brute_optimizer.parameter_grid[710]

print 'actual parameters:', test_parameters
print 'optimizer output: ', global_brute_optimizer(test_signal, parameter_scale_normalization=False)

actual parameters: [  1.55000000e-09   1.58994250e+00  -2.91158829e+00   1.00000000e-10]
optimizer output:  [  1.55000000e-09   1.58994250e+00  -2.91158829e+00   1.00000000e-10]

In this case of course we get the perfect result. However, for actual data this optimizer only returns precomputed brute-force parameter combinations and does not refine them further. It must be given to the Brute2FineOptimizer to get a better result.

Brute Force with refining¶

The Brute2FineOptimizer does not precompute anything, and will do the brute force and refining for every voxel. For this reason it is useful when you want to find specific parameter solutions that are close to some initial parameter guess, but not very efficient when the same global grid can be used for every voxel.

In [6]:

brute2fine_optimizer = brute2fine.Brute2FineOptimizer(mc_model, scheme)

As a simple illustration you can give an initial parameter guess x0 and the optimizer will find a local minimum from that point.

In [7]:

test_signal = global_brute_optimizer.signal_grid[90]
test_parameters = global_brute_optimizer.parameter_grid[90]
x0_vector = global_brute_optimizer.parameter_grid[400]

print 'actual parameters:', test_parameters
print 'start parameters: ', x0_vector
print 'optimizer output: ', brute2fine_optimizer(test_signal, x0_vector) * mc_model.scales_for_optimization

actual parameters: [  2.27500000e-09   1.00507992e+00   3.40817407e-01   1.00000000e-10]
start parameters:  [  1.00000000e-10   8.73750549e-01  -8.94552392e-01   1.00000000e-10]
optimizer output:  [  2.27499914e-09   1.00508076e+00   3.40817944e-01   1.00000000e-10]

Parts of the x0_vector can also be set to None, in which case the optimizer will first do brute-force for those parameters and then do gradient descent to a local minimum.