Applications of Convex Optimization

Rob Knapp, Wolfram Research

ABSTRACT

Convex optimization can be thought of as an extension of linear optimization (linear programming) to nonlinear cases.  Just as the advent of efficient solvers for linear optimization led to and explosion of the applications and uses of linear optimization, the relatively new efficient solvers for convex optimization is opening up a whole new set of applications and uses for convex optimization.  This presentation will describe some of the basic of convex optimization and focus on a variety of examples of applications of convex optimization.

Optimization Classes

Mathematical optimization is typically categorized in 3 separate ways:

Linear/Nonlinear

Constrained/Unconstrained

Local/Global

Taking into consideration convex optimization, a new way to look at the categories is

Convex/General Nonlinear

Constrained/Unconstrained

Local/Global

Actually, for Convex Optimization, the latter two categories are pretty much irrelevant since ...

For a convex optimization problem
                   ({{minimize, f(x)}, {subject to, (Subscript[g, i])(x)≤0}, { , h(x)=0}})
with f(x) and (Subscript[g, i])(x) convex functions and h(x) linear ....

A local optimum is global

The objective may always be considered as linear by transformation (called the “epigraph transformation”) to

({{minimize, t}, { , f(x)≤ t}, {subject to, (Subscript[g, i])(x)≤0}, { , h(x)=0}})

The main reason convex optimization has become more mainstream is the availability of good solvers. (Alas, if you didn’t go, you have missed Nina’s talk on effective solvers, but I’m sure you can still get the notebook.)  These solvers can reliably solve problems that can be expressed in terms of proper convex cones (for information on these go to the reference page for (RefLink[ConicOptimization, paclet:ref/ConicOptimization]) (I need to get to applications at some point!))

Key points that make convex optimization useful for applications:

Optimum is global

Most problems are “solved” once transformed into convex cones

Conic problems have a dual problem that gives sensitivity amongst other things

Applications (at last)

Here is a quick summary of the applications to be demonstrated.

Finance (Portfolio optimization)

Data fitting (L1, cardinality constrained least squares)

Signal processing (Total variation regularization)

Image Processing (Image inpainting)

Geometry (Bounding Ellipsoid, Cone shape)

Graph theory (Max-Cut problem)

Portfolio Optimization

Assuming the mean expected returns and covariances of a set of investments are known, find the mix of investments that balances the return with the variance (risk) of the mix.

{μ,Σ}=;
TableForm[Table[{(IndexedDefault[μ, i]),(IndexedDefault[Σ, i, i])},{i,6}],TableHeadings->{Range[6],{"μ", "σ"}}]

Let (Subscript[w, i]) be the fraction of total investment made in stock i. The return is given by μ.w, where μ is a vector of expected return value of each individual stock:

return=μ.w;

The risk is given by α (w.Σ.w), and α is a risk-aversion parameter:

risk=w.Σ.w;

The objective is to maximize return while minimizing risk for a specified risk-aversion parameter:

objective[α_]:=-(return-α risk);

The percentage of investment w must be greater than 0 and must add to 1:

weightConstraints={w\[VectorGreaterEqual]0,Total[w]==1};

Suppose I fit in with α=1, here is the investment I should make:

QuadraticOptimization[objective[1],weightConstraints,w]

Or if I am more risk-averse

QuadraticOptimization[objective[π],weightConstraints,w]

The solution is fast enough we can easily compute the best mix for a range of risk aversion:

best =Table[Tooltip[{risk,return}/.QuadraticOptimization[
objective[α],weightConstraints,w],ToString[α]],{α,10^Range[-1,1,.1]}];

Here is a plot showing the return of thae best mix along with the returns for 1663 randomly chosen mixes:

randomMixes = Map[((#)/(Total[#]))&, RandomReal[1,{1663, 6}]];
rr = Map[{#.Σ.#,μ.#}&, randomMixes];
Show[ListPlot[best],ListPlot[rr,PlotStyle->Red]]

Data Fitting

Consider finding a cubic polynomial fit for equally spaced points on x∈[-3,3] with given response vector y.

y = {7.,0.,0.,0.,1.,0.,7.};
data = Transpose[{Range[-3.,3.],y}];
basis = {1,x,x^2,x^3};
l2fit = Fit[data,basis,x]

Plot[l2fit,{x,-3,3},Epilog->{PointSize[.025],Point[data]},PlotRange->{-1,8}]

Fit is just using linear algebra to find the least-squares fit.   The design matrix for the basis {1,x,(x^2),(x^3)} with spacing 1 is

A= (({{1, -3, 9, -27}, {1, -2, 4, -8}, {1, -1, 1, -1}, {1, 0, 0, 0}, {1, 1, 1, 1}, {1, 2, 4, 8}, {1, 3, 9, 27}}));

and the standard least-squares fit is

ls = LeastSquares[A,y]

L1 fitting

An issue with least-squares fitting is that deviations are penalized by the square so the fit curve tends to be bent by outliers.  One alternative is to use the L1 norm that penalizes deviations only linearly.  This can be expressed as a convex (actually linear) optimization problem.

l1rules = LinearOptimization[t,Norm[Inactive[Plus][A.c,-y],1]<= t,{c,t}]

Plot[{l2fit,(c /. l1rules).basis},{x,-3,3},Epilog->{PointSize[.025],Point[data]},PlotRange->{-4,8}]

Lp Fitting

In fact, the general Lp norm can be solved by using ConicOptimization (this automatically does a transformation to using a transformation to the “Power Cone”

lptab = Table[( c/. ConicOptimization[t,Norm[Inactive[Plus][A.c,-y],p]<= t,{c,t}]).basis,{p,{1,1.5,2,3,5,8,∞}}];

Plot[Evaluate[lptab],{x,-3,3},Epilog->{PointSize[.025],Point[data]},PlotRange->{-2,8}]

Note: if you are not fond of doing the epigraph transformation yourself,

NMinimize[Norm[Inactive[Plus][A.c,-y],7],c]

automatically uses ConicOptimization

Also you can use Fit directly:

Fit[data,basis,x,NormFunction->7]

Cardinality constrained fitting

Suppose you want to limit the number of basis elements used (i.e. the cardinality or number of nonzero elements of the coefficient vector,  sometimes called the L0 norm).  This can be solved using mixed-integer convex optimization with a binary weight vector.

Note: Paritosh will be giving a presentation at 2PM today specifically on mixed integer convex optimization.

nAllowed = 1;
ccRules = QuadraticOptimization[Norm[Inactive[Plus][A.c,-y]]^2,{Total[w]≤ nAllowed, 0\[VectorLessEqual]w\[VectorLessEqual]1,-bound w\[VectorLessEqual]c\[VectorLessEqual]bound w,bound==100},{c,Element[w,Vectors[4,Integers]]}]

Plot[{l2fit,(c /. ccRules).basis},{x,-3,3},Epilog->{PointSize[.025],Point[data]},PlotRange->{-4,8}]

L1 Regularization

A common thing that is done with convex optimization is to solve a “relaxed” version of a harder problem.  Typically the relaxed version gives the same result most of the time as a harder problem.

Using a regularizing term that is the L1 norm of the coefficients turns out to be a relaxed version of the cardinality constrained problem.

l1reg[λ_] := NArgMin[Norm[Inactive[Plus][A.c,-y]]^2 + λ Norm[c,1],c,Method->"Convex"];

l1reg[1]

An effective way to determine the number of zeros is to use Threshold and SparseArray:

nz[λ_]:= Length[SparseArray[Threshold[l1reg[λ],1.*^-4]]["NonzeroValues"]]

ListPlot[Table[{λ,nz[λ]},{λ,1,6,.1}]]

l1reg[6]

Basis Pursuit

This example is actually sneaking into the territory of signal processing, but it is just a bigger example of the use of L1 regularization.

Given a set of basis functions (Subscript[ϕ, i])(t), a signal x(t) may be approximated by the linear combination x(t)≃(∑_(j=1)^n)(Subscript[α, j])(Subscript[ϕ, j])(t). Typically, the larger the basis, the better a signal can be approximated. However, not all signals need all of the basis functions for a reasonable approximation. The idea of basis pursuit is to choose a small subset of a larger set, or dictionary, of basis functions that may be overdetermined.

Use a basis of Gabor functions (e^(-(((((t-(Subscript[t, 0]))/σ)))^2)))cos(ω (t-(Subscript[t, 0]))) and (e^(-(((((t-(Subscript[t, 0]))/σ)))^2)))sin(ω (t-(Subscript[t, 0]))) for an example of a large dictionary of basis functions.

cosGabor[t0_,ω_, σ_, t_] := Exp[-((t-t0)/σ)^2] Cos[ω (t-t0)];
sinGabor[t0_,ω_, σ_, t_] := Exp[-((t-t0)/σ)^2] Sin[ω (t-t0)];

Here is a plot of showing a sample of Gabor functions.

Plot[{cosGabor[.25,10,.05,t],sinGabor[.25,10,.05,t],cosGabor[.5,50,.05,t],sinGabor[.5,50,.05,t],cosGabor[.75,150,.05,t],sinGabor[.75,150,.05,t]},{t,0,1},PlotRange->All, AspectRatio->.3,ImageSize->Large]

On the interval from 0 to 1 with a time sample interval δt, a set of basis functions with (Subscript[t, 0])=0,δt,2δt,…,1, ω=0, 5, 10,…,150 and a fixed σ=1/20 are sufficient to exactly represent a wide variety of signals.

σ = 0.05;
nt = 500;
ts =Subdivide[0.,1.,nt];
δt = ts[[1]];

Length[basis = Flatten[Table[{cosGabor[t0,0.,σ,t],Table[{cosGabor[t0,ω,σ,t],sinGabor[t0,ω,σ,t]},{ω,5,150,5}]},{t0,ts}]]]

The basis matrix can be constructed using DesignMatrix.

Dimensions[designMatrix = DesignMatrix[Thread[{ts,0.}],basis, t,IncludeConstantBasis->False]]

This matrix is very large, but with a small loss in representation accuracy, it is possible to use the fact the the Gabor functions decay quickly from (Subscript[t, 0]) to make it sparse by using a threshold.

basisMatrix = SparseArray[Threshold[designMatrix, 0.001]]

Now consider a signal.

signal[t_] := Exp[-9(t-.5)^2](1 + .5 Sin[13t])Sin[27 Sin[7t]];
Plot[signal[t],{t,0,1}]

The signal can be discretized by just evaluating at all the time samples.

discreteSignal = signal[ts];

Now use L1 regularization to find a sparse fit.

l1Fit = Fit[{basisMatrix, discreteSignal},FitRegularization->{"L1", 1.}]

The sparse representation uses only 53 out of 30561 basis elements.  Note that the solution to this problem gives the coefficients that minimize the sum of the L2 norm of the residuals and the L1 term.  What we want is to find the coefficients that minimize the L2 norm of the residuals for the chosen basis elements.

chosen = Flatten[l1Fit["NonzeroPositions"]];
dc = basisMatrix[[All,chosen]];
chosenL2 = Fit[{dc, discreteSignal}];
sparseFitFun[t_] = chosenL2.basis[[chosen]];

Show the error.

{Plot[{signal[t],sparseFitFun[t]},{t,0,1}], Plot[signal[t]-sparseFitFun[t],{t,0,1},PlotRange->All]}

There is a tradeoff between the error and the number of basis functions used.  Incresing λ decreases the number of basis functions and increases the error.

l1Fit = Fit[{basisMatrix, discreteSignal},FitRegularization->{"L1", 1.}]

The sparse representation uses only 53 out of 30561 basis elements.  Note that the solution to this problem gives the coefficients that minimize the sum of the L2 norm of the residuals and the L1 term.  What we want is to find the coefficients that minimize the L2 norm of the residuals for the chosen basis elements.

chosen = Flatten[l1Fit["NonzeroPositions"]];
dc = basisMatrix[[All,chosen]];
chosenL2 = Fit[{dc, discreteSignal}];
sparseFitFun[t_] = chosenL2.basis[[chosen]];

Show the error.

{Plot[{signal[t],sparseFitFun[t]},{t,0,1}], Plot[signal[t]-sparseFitFun[t],{t,0,1},PlotRange->All]}

There is a tradeoff between the error and the number of basis functions used.  Incresing λ decreases the number of basis functions and increases the error.

l1Fit = Fit[{basisMatrix, discreteSignal},FitRegularization->{"L1", 1.}]

The sparse representation uses only 53 out of 30561 basis elements.  Note that the solution to this problem gives the coefficients that minimize the sum of the L2 norm of the residuals and the L1 term.  What we want is to find the coefficients that minimize the L2 norm of the residuals for the chosen basis elements.

chosen = Flatten[l1Fit["NonzeroPositions"]];
dc = basisMatrix[[All,chosen]];
chosenL2 = Fit[{dc, discreteSignal}];
sparseFitFun[t_] = chosenL2.basis[[chosen]];

Show the error.

{Plot[{signal[t],sparseFitFun[t]},{t,0,1}], Plot[signal[t]-sparseFitFun[t],{t,0,1},PlotRange->All]}

There is a tradeoff between the error and the number of basis functions used.  Incresing λ decreases the number of basis functions and increases the error.

l1Fit = Fit[{basisMatrix, discreteSignal},FitRegularization->{"L1", 10.}]

chosen = Flatten[l1Fit["NonzeroPositions"]];
dc = basisMatrix[[All,chosen]];
chosenL2 = Fit[{dc, discreteSignal}];
sparseFitFun[t_] = chosenL2.basis[[chosen]];

{Plot[{signal[t],sparseFitFun[t]},{t,0,1}], Plot[signal[t]-sparseFitFun[t],{t,0,1},PlotRange->All]}

Signal Processing

There are many ways convex optimization can be used for signal processing.  I’ll just show one simple example.

Total Variation Regularization

For a signal (Subscript[x, i]), this is done by finding values (Subscript[s, i]) that minimize ‖s-x‖+λ(Subscript[‖δs‖, 1]), where δs is the vector with elements (Subscript[s, i])-(Subscript[s, i-1]),i = 2,…,n and λ is a parameter that is used to trade off between noise reduction and proximity to the signal x.

Consider the Dow Jones Industrial Average stock index during the 2008 financial crisis. There is a lot of day-to-day volatility that over a long term can effectively be considered noise.

DateListPlot[index=TimeSeries[FinancialData["^DJI",{"Jan. 1, 2008","Dec. 31 2010"}]]]

The signal is the sequence of values of the stock index.

values = QuantityMagnitude[index["Values"]];
n = Length[values];

Find the smoothed index for λ=1. Method->"Convex" is used, since other methods would be very slow if they succeeded at all for this many variables.

s[1]=Quiet[NArgMin[{Norm[s-x] + Norm[Differences[s],1],x==values},Element[s,Vectors[n]],Method->"Convex"]];

Compare the smoothed index with the actual index.

DateListPlot[{index, TimeSeries[Transpose[{index["Dates"],QuantityArray[s[1], "USDollars"]}]]},PlotStyle->{Automatic, Thickness[0.01]},ImageSize->Large]

It is interesting to look at the smoothed index for some of the different values of λ.

Grid[Map[ListLinePlot[s[#],Axes->False,PlotLabel->"λ="<>ToString[#,TraditionalForm]]&,Partition[2^Range[-5,3],3],{2}]]

Image Processing

Image Inpainting

The idea here is to try to reconstruct a corrupted image.  Suppose that you have only a subset of pixels of an image.  Inpainting minimizes the norms of the variation between neighboring pixels to make an informed guess.

Here we manually corrupt a good image for demonstration.

dir = NotebookDirectory[];
ImageDimensions[image = ColorConvert[Import[FileNameJoin[{dir,"Mystery.jpg"}]],"Grayscale"]]

data = ImageData[image];
{m, n} = Dimensions[data];
BlockRandom[
mask = SparseArray[RandomChoice[{0.6,0.4}->{0,1},Dimensions[data]]],RandomSeeding->123];
corruptedData=mask*data;
knownPositions=mask["NonzeroPositions"];
Image[corruptedData]

The convex optimization problem is 
           ({{minimize, (∑_(i=1)^(m-1))(∑_(j=1)^(n-1))(Subscript[(Norm[{(Subscript[u, i+1,j])-(Subscript[u, i,j]),(Subscript[u, i,j+1])-(Subscript[u, i,j])}]), 2])}, {subject to , (Subscript[u, i,j])=(Subscript[data, i,j]) ∀{i,j}∈knownPositions}})

This can be done with NMinimize but it is more efficient (this is something we are working on, but has to do with the way symbols are processed in WL) using ConicOptimization

normConstraints = Table[Inactive[Norm][{Indexed[u,{i+1,j}]-Indexed[u,{i,j}], Indexed[u,{i,j+1}]-Indexed[u,{i,j}]}]≤ Indexed[t,{i,j}],{i,m-1},{j,n-1}];

equalityConstraints = Map[Function[Indexed[u,#]==Extract[data, #]], knownPositions];

AbsoluteTiming[rule = ConicOptimization[Total[t,2],{normConstraints, equalityConstraints},{Element[u, Matrices[{m,n},Reals]], Element[t,Matrices[{m-1,n-1},Reals]]},Tolerance->1.*^-3];]

Image[u /. rule]

image

Geometry

Many geometrical problems can be mapped to convex optimization.  Here I’ll just show two examples

Minimal Bounding Ellipsoid

Optimal Ice Cream Cone

Minimal Bounding Ellipsoid

Given a set of m points (Subscript[x, i]) in n dimensions, find the Löwner–John ellipsoid that is the smallest-volume ellipsoid enclosing all of the points.

This example demonstrates that for problems expressed in terms of highly nonlinear convex objective functions and constraints, NMinimize can do transformations that make the problem easily solved by convex optimization.

An ellipsoid in (ℝ^n) can be described by the set of points x such that ‖A.x+b‖≤1, where A is a positive semidefinite n×n matrix and b is an n-vector. The volume of the ellipsoid is proportional to det((A^(-1))). log(det((A^(-1))))=-log(det(A)) is convex so finding the solution of

({{minimize,  -log(det(A))}, {subject to,             ({{‖A.(Subscript[x, i])+b‖≤ 1, i=1… m}, {A(Underscript[\[VectorGreaterEqual], (SubPlus[S])])0}})}})

The example shown here is done for a random set of points in three dimensions so that it can be visualized, but the solution algorithm can handle larger dimensions n.

n = 3;
points = BlockRandom[RandomReal[{-1,1},{47,3}],RandomSeeding->47];

Find the solution using NMinimize.

{minval, rules}= NMinimize[{-Log[Det[a]],{a(Underscript[\[VectorGreaterEqual], (SemidefiniteConeString)])0, Table[Norm[a.xi+b]≤1,{xi,points}]}}, {a, b∈Vectors[n]}]

In three dimensions, the ellipsoid can be described in terms of the Ellipsoid graphics primitive.

ainv = Inverse[a /. rules];
center = -ainv.(b /. rules);
Σ = ainv.ainv;
RegionMeasure[ellipsoidLJ = Ellipsoid[center,Σ]]

ellipsoidFast = BoundingRegion[points, "FastEllipsoid"];
RegionMeasure[ellipsoidFast]

Graphics3D[{PointSize[0.025], Point[points], Opacity[.25],{Green, ellipsoidLJ},{LightBlue,ellipsoidFast} },Axes->True]

Optimal Ice Cream Cone

Design an ice cream cone for a spherical scoop of ice cream of volume V. Assume the radii of the scoop and cone are the same.  Minimize the surface area the cone so that it cone would hold the entire scoop when melted:

Minimize the surface area of and ice cream cone so that the cone would hold at least an entire spherical scoop of ice cream of volume V the same radius when melted:

coneArea =π r (r+(Sqrt[(h^2)+(r^2)]));

coneVolume = (π/3) h (r^(2 ));

scoopVolume = (4/3)π (r^3);

Now a new function ... GeometricOptimization.  This takes problems expressed in terms of sums of products of powers of positive terms (not necessarily convex) and transforms to a convex form (using a logarithmic transformation) to a convex form.

GeometricOptimization[coneArea, {scoopVolume ≤ coneVolume, scoopVolume == V, V == 1}, {h, r}]

Note that it is NOT called GeometricOptimization because it is only for geometric problems.

Graph Theory

MaxCut problem

The max-cut problem determines a subset S of the vertices V of a graph, for which the sum of the weights (Subscript[w, i j]) of the edges that cross from S to its complement V \ S is maximized.

The max-cut problem is NP-complete, but a relaxation of the max-cut problem can be solved efficiently using semidefinite optimization.

Suppose that V has n vertices so that they can be described by an index 1≤k≤n. Let x be a vector with components (Subscript[x, k])=1 for k∈S and (Subscript[x, k])=-1 for k∈V \ S so that 1-(Subscript[x, i]) (Subscript[x, j]) is nonzero (=2) only if i∈S and j∈V \ S. Thus the max-cut is found by maximizing (1/4)(∑_(i,j=1)^n)(Subscript[w, i j ])(1-(Subscript[x, i])(Subscript[x, j])), where (Subscript[x, k])∈{-1, 1}, k=1,…, n. The objective can be rewritten as:

(1/4)(∑_(i,j=1)^n)((Subscript[x, i])(Subscript[w, i j ])(Subscript[x, i])-(Subscript[x, i])(Subscript[w, i j ])(Subscript[x, j]))=(1/4)xᵀ·(W·(1, 1, …,1) - W)·x =(1/4) xᵀ·L·x =(1/4)Tr ( L · (x·xᵀ))

... where L is the Laplacian matrix of the graph and W is the weighted adjacency matrix.

laplacianMatrix[G_?GraphQ]:=(DiagonalMatrix[Total[#,{1}]]-#)&[WeightedAdjacencyMatrix[G]]

For smaller cases, the max-cut problem can be solved exactly, but this is impractical for larger graphs since in general the problem has NP-complete complexity.

maxCutExact[G_?GraphQ]:=Module[{L=laplacianMatrix[G],n,x},
n=Length[First[L]];
x=Array["var",n];
MaxValue[{1/4x.L.x,x^2==ConstantArray[1,n]},x,Integers]
]

The problem minimizes (1/4) Tr[L·Y], where Y=x·(Transpose[x]) is a symmetric rank-1 positive semidefinite matrix, with (Subscript[Y, i i])=1 for each i, equivalent to Tr[(Subscript[E, i ])·Y]=1, where (Subscript[E, i]) is the matrix with 1 at the i(^th) diagonal position and 0 everywhere else.

e[i_,n_] := SparseArray[{i,i}->1.,{n,n},0.];

To make the solution practical, solve a relaxed problem where the rank-1 condition is eliminated so that it only requires that Y be positive semidefinite.  The (relaxed) convex problem is expressed as
            ({{maximize, -Tr[L.y]}, {subject to , Tr[(Subscript[e, i]).y]=1,i=1,…,n}, { , y\[VectorGreaterEqual]0}})

This is the dual problem for the semidefinite problem
                      ({{minimize, {1,…,1}.v}, {subject to, -(1/4)L+(∑_(i=1)^n)(Subscript[e, i])(Subscript[v, i])\[VectorGreaterEqual]0}})

For the solution Y of the relaxed problem, a cut is constructed by randomized rounding: decompose Y=U (Transpose[U]), let ρ be a uniformly distributed random vector of the unit norm and let (Subscript[x, k])=1⟺ρ·(Subscript[U, k])>0. For demonstration, a function is defined that shows the relaxed value, the rounded value and the graph, with the vertices in S shown in red.

maxCutRelaxed[G_?GraphQ, trials_:10]:=Module[{L=laplacianMatrix[G],Y,U,n,sdpmax,x},
n=Length[First[L]];{Y,sdpmax}=SemidefiniteOptimization[ConstantArray[1,n], Join[{-L/4},Table[e[i,n], {i,n}]],{"DualMaximizer","DualMaximumValue"}];
U=CholeskyDecomposition[Y];
x=First[MaximalBy[Table[Sign[RandomPoint[Sphere[n]].U],trials],1/4#.L.#&]];
{sdpmax,1/4x.L.x, HighlightGraph[G,Pick[Range[n],x,1]]}
]

Lets test it on an actual graph

pG=PetersenGraph[13,5,VertexSize->Medium]

maxCutRelaxed[pG]

Summary

I have shown a very few applications of convex optimization that at least begins to show the breadth of what can be done taking advantage of reliable solvers for convex problems.

There are many more applications shown in the reference pages for the Wolfram Language optimization functions:

(RefLink[LinearOptimization, paclet:ref/LinearOptimization])

(RefLink[LinearFractionalOptimization, paclet:ref/LinearFractionalOptimization])

(RefLink[QuadraticOptimization, paclet:ref/QuadraticOptimization])

(RefLink[SecondOrderConeOptimization, paclet:ref/SecondOrderConeOptimization])

(RefLink[SemidefiniteOptimization, paclet:ref/SemidefiniteOptimization])

(RefLink[ConicOptimization, paclet:ref/ConicOptimization])

(RefLink[GeometricOptimization, paclet:ref/GeometricOptimization])