MATLAB Examples

BertiniLab: User Guide

The BertiniLab package is designed to convert numerical MATLAB® functions into text files that can be used as input to the numerical algebraic geometry package Bertini. This package can find all the roots of systems of polynomial equations. The package includes two classes: polysym and BertiniLab.

The class polysym implements polynomial calculations in a symbolic form. It is similar to the sym class in the MATLAB symbolic toolbox, but tailored to the requirements of Bertini. If the inputs to a function are defined as polysym objects, the output will also be a polysym object.

Given a function and arguments, definitions and configurations defined as polysym objects, BertiniLab can be used to create the input file for Bertini, run Bertini, and collect the outputs.

This guide concentrates on how to use BertiniLab as an interface with Bertini. For how Bertini is used, see Numerically solving polynomial systems with Bertini, by Daniel J. Bates, Jonathan D. Haunstein, Andrew J. Sommese and Charles W. Wampler (SIAM 2013). The examples from the book are also implemented in BertiniLab.

Contents

Creating a Bertini object

An empty Bertini object is created by the call

poly_system = BertiniLab;
disp(BertiniLab)
  BertiniLab with properties:

                 config: [0x0 struct]
           function_def: []
               variable: []
         variable_group: []
     hom_variable_group: []
           pathvariable: []
                 random: []
            random_real: []
     definedSubfunction: []
              parameter: [0x2 polysym]
               constant: [0x2 polysym]
            subfunction: [0x2 polysym]
        starting_points: [0x0 struct]
       final_parameters: []
          member_points: []
    interactive_choices: {}
      Bertini_io_folder: ''
          view_progress: 0
          solve_summary: ''
          function_name: []

Most of the properties correspond to Bertini declarations. In particular, variable_group is the most common property for assigning variables. These can be assigned in a statement like

poly_system.variable_group = 'x';

Similarly, the equation(s) to be solved are assigned to the property function_def:

poly_system.function_def = 'x^2-1';

Equivalently, these properties can be defined in the initial call to BertiniLab using name/value pairs:

poly_system = BertiniLab('variable_group','x','function_def','x^2-1'); %#ok<*NASGU>

Defining the variables

Bertini solves a system of polynomial equations, so it needs a list of variables and the equations in those variables.

Depending on the problem, the following properties can be used for defining variables:

As an illustration, we will solve a generalized eigenvalue problem

$$\mu A v = \lambda B v. $$

If the matrices are 2x2, the variables are $\mu$, $\lambda$, $v_1$ and $v_2$. For BertiniLab, these can be defined as a cell array of strings:

vars = {'mu0','lambda','v1','v2'};

(There is a function mu in a MATLAB toolbox, so to avoid unexpected behavior we will call the variable mu0).

Alternatively, the variables can be defined as polysym objects:

polysyms mu0 lambda v1 v2
vars = [mu0 lambda v1 v2];

Using polysym, the components of v can be defined in vector form:

v = polysym('v',[1 2]);
vars = [mu0 lambda v];

In this example, the equations are homogeneous, so the variables are assigned to the property hom_variable_group:

poly_system = BertiniLab('hom_variable_group',vars);

It is better to use multihomogeneous homotopy to solve this problem. This involves grouping variables in a way that reduces the number of paths that Bertini must track. To group variables represented by polysym objects, we group them into vectors and put the result in a cell array:

poly_system = BertiniLab('hom_variable_group',{[mu0 lambda],v});

(The two variables in v are already in a vector.)

BertiniLab also accepts variable precision and symbolic objects.

Defining the equations

Equations are defined in the property function_def. We will look at a specific example of the matrices:

$$A = \left[\begin{array}{cc} 1 & 2 \\ 2 & 4 \end{array}\right]$$

$$B = \left[\begin{array}{cc} 4 & -2 \\ -2 & 1 \end{array}\right].$$

Written out in full, the resulting equations are

$$ \mu \left(v_1 + 2 v_2\right) - \lambda \left(4 v_1 - 2 v_2\right) = 0 $$

$$ \mu \left(2 v_1 + 4 v_2\right) = \lambda \left(-2 v_1 + v_2\right)= 0. $$

There are three basic ways to define the equations for BertiniLab to solve:

  1. Write out the equations explicitly
  2. Use a pre-defined function
  3. Create a subclass

Write out the equations explicitly

One way to define the equations is as a cell array of strings:

eqns = {'mu0*(v1+2*v2) - lambda*(4*v1-2*v2)', 'mu0*(2*v1+4*v2) - lambda*(-2*v1+v2)'};

If the variables are polysym objects, the functions can be defined using matrix algebra:

A = [1 2; 2 4]; B = [4 -2; -2 1];
eqns = mu0*A*v.' - lambda*B*v.';

Using a pre-defined function

The equations can also be defined using a function call. This function can be pre-defined. For example, suppose we have the function

f = @(X,Y,m,l,w) (m*X-l*Y)*w.';

We can now define the variables and produce expressions for the equations.

polysyms mu0 lambda
A = [1 2; 2 4]; B = [4 -2; -2 1];

A and B will be determined by the physical problem. Here we'll use the same A and B, but formulate the problem so that larger matrices could be substituted.

N = length(A);
v = polysym('v',[1 N]);
eqns = f(A,B,mu0,lambda,v);

The Bertini input file expects each equation to have a name. By default, BertiniLab generates the names. If the user wishes to provide the names, they can use the property function_name. For example, to call the equations f1, f2, etc., use the following:

poly_system.function_name = polysym('f',[1 N]);

Solving BertiniLab

We will now combine the variable and function definitions and solve for B

vars = {[mu0 lambda],v};
poly_system = BertiniLab('function_def',eqns,'hom_variable_group',vars);
poly_system = solve(poly_system);

The above code creates an input file for Bertini with the following content:

type('./input')
CONFIG
END;

INPUT
hom_variable_group mu0, lambda;
hom_variable_group v1, v2;
function eqn1, eqn2;
eqn1=(mu0-lambda*4)*v1+(mu0*2-lambda*(-2))*v2;
eqn2=(mu0*2-lambda*(-2))*v1+(mu0*4-lambda)*v2;
END;

Of course, for this example it would be just as easy to create the file by hand; but in creating it, BertiniLab also stores some information that makes it easy to import the solutions into MATLAB and analyze them.

After creating the input file (and any other necessary data files), the solve method runs Bertini, which creates several output files.

Note that you need to assign the output of the solve to poly_system. This saves information that is needed to access and interpret the output of Bertini. In addition, the property solve_summary contains the diagnostic output from Bertini:

disp(poly_system.solve_summary)
    Bertini(TM) v1.4
   (October 23, 2013)

 D.J. Bates, J.D. Hauenstein,
 A.J. Sommese, C.W. Wampler

(using GMP v5.1.3, MPFR v3.1.2)



NOTE: You have requested to use adaptive path tracking.  Please make sure that you have
setup the following tolerances appropriately:
CoeffBound: 5.400840000000e+00, DegreeBound: 2.000000000000e+00
AMPSafetyDigits1: 1, AMPSafetyDigits2: 1, AMPMaxPrec: 1024

Tracking path 0 of 2

Solution Summary

NOTE: nonsingular vs singular is based on condition number and identical endpoints

		| Number of real solns	|  Number of non-real solns	|  Total
------------------------------------------------------------------------------------------
Non-singular	|	2		|		0		|   2
Singular	|	0		|		0		|   0
------------------------------------------------------------------------------------------
Total		|	2		|		0		|   2


Multiplicity Summary

  Multiplicity	|  Number of real solns	|  Number of non-real solns
------------------------------------------------------------------------------------------
	1	|	2		|	0


------------------------------------------------------------------------------------------
The following files may be of interest to you:

main_data:             A human-readable version of the solutions - main output file.
raw_solutions:         A list of the solutions with the corresponding path numbers.
raw_data:              Similar to the previous, but with more information about the solutions.
real_solutions:        A list of all real solutions.
nonsingular_solutions: A list of all nonsingular solutions.
singular_solutions:    A list of all singular solutions.
------------------------------------------------------------------------------------------

Paths Tracked: 2

Retrieving the solutions

The method read_raw_data reads the diagnostic information for the solutions. One solution is displayed below.

diagnostics = poly_system.read_raw_data;
disp(diagnostics(1))
                                 path_number: 0
                               max_precision: 52
                                 coordinates: [1x2 polysym]
                           function_residual: 9.9301e-16
                            condition_number: 15.8136
                       final_Newton_residual: 3.2418e-16
                          final_pathvariable: 3.9063e-04
                                    lastdiff: 6.2420e-15
    pathvariable_of_first_precision_increase: 0
                        endgame_cycle_number: 1
                                   path_flag: 1

The method read_solutions can be used to read most other data files:

sols = poly_system.read_solutions('nonsingular_solutions');
disp(sols)
  1.000000000000000e+00+(0.000000000000000e+00)*I -8.127142098938020e-18+(-1.031326620820564e-17)*I
                                                                                                   
-2.538891751908466e-17+(-9.838741889091606e-19)*I   1.000000000000000e+00+(0.000000000000000e+00)*I
                                                                                                   
  1.000000000000000e+00+(0.000000000000000e+00)*I   5.000000000000001e-01+(0.000000000000000e+00)*I
                                                                                                   
-4.999999999999999e-01+(-5.551115123125783e-17)*I   1.000000000000000e+00+(0.000000000000000e+00)*I

The solutions are saved as polysym arrays, so they retain the full accuracy of the Bertini solutions. The class polysym has the method real_imag to extract the real and imaginary parts of the solution:

[re,im] = sols.real_imag;
disp(re)
 1.000000000000000e+00 -8.127142098938020e-18
                                             
-2.538891751908466e-17  1.000000000000000e+00
                                             
 1.000000000000000e+00  5.000000000000001e-01
                                             
-4.999999999999999e-01  1.000000000000000e+00

Solutions can also be converted to double precision, variable precision or symbolic objects.

disp(re.double)
    1.0000   -0.0000
   -0.0000    1.0000
    1.0000    0.5000
   -0.5000    1.0000

disp(re.vpa)
[                                 1.0, -0.00000000000000000812714209893802]
[ -0.00000000000000002538891751908466,                                 1.0]
[                                 1.0,                  0.5000000000000001]
[                 -0.4999999999999999,                                 1.0]
 

The data can also be matched with the variables using match_solutions. This puts the data into a structure with a field for each variable, named after that variable. Vector and matrix variables are collected in arrays of appropriate size.

sols = poly_system.match_solutions('nonsingular_solutions');
disp(sols)
disp(sols.v)
    lambda: [2x1 polysym]
       mu0: [2x1 polysym]
         v: [2x2 polysym]

  1.000000000000000e+00+(0.000000000000000e+00)*I 5.000000000000001e-01+(0.000000000000000e+00)*I
                                                                                                 
-4.999999999999999e-01+(-5.551115123125783e-17)*I 1.000000000000000e+00+(0.000000000000000e+00)*I

This structure can be converted back to polysym arrays using unpack_solutions.

sols = poly_system.unpack_solutions(sols);
disp(sols)
  1.000000000000000e+00+(0.000000000000000e+00)*I -8.127142098938020e-18+(-1.031326620820564e-17)*I
                                                                                                   
-2.538891751908466e-17+(-9.838741889091606e-19)*I   1.000000000000000e+00+(0.000000000000000e+00)*I
                                                                                                   
  1.000000000000000e+00+(0.000000000000000e+00)*I   5.000000000000001e-01+(0.000000000000000e+00)*I
                                                                                                   
-4.999999999999999e-01+(-5.551115123125783e-17)*I   1.000000000000000e+00+(0.000000000000000e+00)*I

This format is needed for the data in some applications where Bertini reads additional data files.

Using a subclass

BertiniLab is a Matlab class, so the user can create subclasses to modify its behavior or apply it to particular kinds of problem. The following code (which can be found in the file generalized_eigenvalue_problem.m) constructs a subclass of BertiniLab that takes just two arguments, the matrices A and B, finds the generalized eigenvalues mu and lambda, and saves them as properties of the subclass. Note that the method solve is overloaded for this purpose.

classdef generalized_eigenvalue_problem < BertiniLab
   properties
      mu = [];
      lambda = [];
   end
   methods
      function obj = generalized_eigenvalue_problem(A,B)
          m = polysym('mu');
          l = polysym('lambda');
          N = length(A);
          v = polysym('v',[1 N]);
          obj = obj@BertiniLab('function_def',(m*A-l*B)*v.','hom_variable_group',{[m l],v});
          obj.mu = m;
          obj.lambda = l;
      end
      function obj = solve(obj)
          obj = solve@BertiniLab(obj);
          mu = obj.mu; lambda = obj.lambda;
          sstruct = obj.match_solutions('nonsingular_solutions','mu','lambda');
          obj.mu = sstruct.mu;
          obj.lambda = sstruct.lambda;
      end
   end
end

Here is the code to run it:

A = [1 2; 2 4]; B = [4 -2; -2 1];

poly_system = generalized_eigenvalue_problem(A,B);
poly_system = poly_system.solve;
disp([poly_system.mu0 poly_system.lambda])
1.000000000000000e+00+(0.000000000000000e+00)*I -3.305883366651087e-17+(-1.283405622485310e-17)*I
                                                                                                 
1.216673641681777e-17+(1.640181098080491e-17)*I   1.000000000000000e+00+(0.000000000000000e+00)*I

Shortcuts and utilities

Some simple shortcuts are available for using BertiniLab.

polysolve has a syntax similar to that of the Symbolic Math ToolBox function solve. For example, to solve the system of equations

$$ x^2-1 = 0 $$

$$ y^2-4 =0, $$

for which the solution set consists of the four points

$$ (\pm 1,\pm 2), $$

the following command can be used:

sols = polysolve({'x^2-1'; 'y^2-4'});
disp(sols)
    x: [4x1 polysym]
    y: [4x1 polysym]

The solutions are

disp([sols.x sols.y])
 1.000000000000000e+00+(0.000000000000000e+00)*I  2.000000000000001e+00+(3.330669073875470e-16)*I
                                                                                                 
1.000000000000000e+00+(-3.469446951953614e-18)*I -2.000000000000000e+00+(6.938893903907228e-18)*I
                                                                                                 
-1.000000000000000e+00+(0.000000000000000e+00)*I  2.000000000000000e+00+(2.220446049250313e-16)*I
                                                                                                 
-1.000000000000000e+00+(0.000000000000000e+00)*I -2.000000000000000e+00+(0.000000000000000e+00)*I

The command findRealRoots has a similar syntax to polysolve but returns the real roots. For example,

sols = findRealRoots({'x^2-1'; 'y^2-4'});
disp([sols.x sols.y])
    1.0000    2.0000
    1.0000   -2.0000
   -1.0000    2.0000
   -1.0000   -2.0000

The command BertiniClean removes all the Bertini-generated files from a folder.

Configuring BertiniLab

Choosing an I/O folder:

BertiniLab creates input files for Bertini and reads output files. By default, all the files go in the current folder (which can be identified using pwd). However, the user can change this using the property Bertini_io_folder. For example, the following commands provide the full path name for a folder into which all the files will go.

poly_system = BertiniLab;
poly_system.Bertini_io_folder = '/Users/my_name/matlab/Bertini_IO';

A relative path name can be provided instead. The following will put all the files in a subfolder Bertini_IO of the current folder (creating Bertini_IO if necessary).

poly_system.Bertini_io_folder = './Bertini_IO';

The above changes only last as long as the BertiniLab object exists. If the user wants to use the same IO folder for all BertiniLab runs, they can edit the line Bertini_io_folder = ''; in the file BertiniLab.m.

Monitoring a Bertini run:

Bertini reports stages in the progress of a run (for example, how many curves have been traced). By default, this information is stored in the variable solve_summary. However, if the sets the property view_progress to true, this information will be echoed to the command window, so the user can monitor the progress of the calculation:

poly_system = BertiniLab;
poly_system.view_progress = true;