diff --git a/@msspoly/diff.m b/@msspoly/diff.m index eeb5c25..dbe4090 100644 --- a/@msspoly/diff.m +++ b/@msspoly/diff.m @@ -35,25 +35,32 @@ p = p(:); end -% Find rows with xn(i) in them -match = msspoly.match_list(xn,p.var); -msk = match ~= 0; -[i,j] = ind2sub(size(p.var),find(msk)); +[q,err]=double(p); -% match(msk) -- which column to place result in -% i -- which entries to draw var/pow/coeff from -% j -- which power to decrement / multiply by coeff - -var = p.var(i,:); -coeff = p.coeff(i).*p.pow(msk); -pow = p.pow(i,:); -ind = sub2ind(size(pow),(1:length(i))',j); - -pow(ind) = pow(ind) - 1; - -J = msspoly([size(p,1) length(xn)],... - [ p.sub(i,1) match(msk) ],... - var,pow,coeff); +if ~err + J = zeros(size(p,1),size(x,1)); + +else + % Find rows with xn(i) in them + match = msspoly.match_list(xn,p.var); + msk = match ~= 0; + [i,j] = ind2sub(size(p.var),find(msk)); + + % match(msk) -- which column to place result in + % i -- which entries to draw var/pow/coeff from + % j -- which power to decrement / multiply by coeff + + var = p.var(i,:); + coeff = p.coeff(i).*p.pow(msk); + pow = p.pow(i,:); + ind = sub2ind(size(pow),(1:length(i))',j); + + pow(ind) = pow(ind) - 1; + + J = msspoly([size(p,1) length(xn)],... + [ p.sub(i,1) match(msk) ],... + var,pow,coeff); +end if nargin == 3 J = reshape(J*d,szp); diff --git a/@spotsqlsol/spotsqlsol.m b/@spotsqlsol/spotsqlsol.m deleted file mode 100644 index 9770363..0000000 --- a/@spotsqlsol/spotsqlsol.m +++ /dev/null @@ -1,51 +0,0 @@ -classdef spotsqlsol - properties - program = []; - primalSolution = []; - dualSolution = []; - dualSlack = []; - objective = []; - dualObjective = []; - info = []; - end - - methods - function f = primalInfeasible(sol) - f = sol.info.pinf; - end - function f = dualInfeasible(sol) - f = sol.info.dinf; - end - - function sol = spotsqlsol(prog,info,obj,dobj,psol,dsol,dslack) - sol.program = prog; - sol.info = info; - sol.primalSolution = psol; - sol.dualSolution = dsol; - sol.dualSlack = dslack; - sol.objective = obj; - sol.dualObjective = dobj; - end - - function ev = eval(sol,exp) - if sol.primalInfeasible - error('Primal Infeasible, cannot evalutate primal variables.'); - end - ev = subs(exp,[sol.program.variables],... - [sol.primalSolution]); - end - function ev = dualEval(sol,exp) - if sol.dualInfeasible - error('Dual Infeasible, cannot evalutate dual variables.'); - end - ev = subs(exp,[sol.program.dualVariables],... - [sol.dualSolution]); - end - function ev = dualSlackEval(sol,exp) - if sol.dualInfeasible - error('Dual Infeasible, cannot evaluate dual variables.'); - end - ev = subs(exp,sol.program.variables,sol.dualSlack); - end - end -end \ No newline at end of file diff --git a/README.txt b/README.txt deleted file mode 100644 index c93ca33..0000000 --- a/README.txt +++ /dev/null @@ -1,5 +0,0 @@ -This is a Subversion repository; use the 'svnadmin' tool to examine -it. Do not add, delete, or modify files here unless you know how -to avoid corrupting the repository. - -Visit http://subversion.tigris.org/ for more information. diff --git a/bin/.gitignore b/bin/.gitignore new file mode 100644 index 0000000..7b58de6 --- /dev/null +++ b/bin/.gitignore @@ -0,0 +1,4 @@ +*.mexmaci64 +*.mexmaci +*.mexa64 +*.mexglx diff --git a/changelog.txt b/changelog.txt deleted file mode 100644 index 19a095a..0000000 --- a/changelog.txt +++ /dev/null @@ -1,73 +0,0 @@ -SPOT change log --------------- - -22.03.10 -Initial SVN commit. - -28.02.10 - -a better working ltid_pass.m -ltid_vwqt2ph.m - -25.02.10 - -mss_m2bs.m to complement mss_v2s.m -ltid_pass.m to replace ltid_passive.m - -11.01.10 - -introduced spot_chk.m - -05.10.09 - -fixed a bug in msspoly/subsasgn.m -new function msspoly/msubs.m for vectorized sampling of msspoly -simulation and piecewise polynomial function handling in the new pot/sim directory -updated pot_install.m - -15.07.09 - -introduced accurate algorithm option in ltid_vw2ab.m and ltid_vw2ab_psd.m - -07.07.09 - -ltid functions re-written and re-named to follow a uniform convention - -01.07.09 - -passive MIMO LTI sysid code, summarized in ltid_passive.m is added - -24.06.09 - -fixed a bug in ltid_abpsd.m - -06.05.09 - -added the option of registering -multiple SeDuMi-style decision variable -cones (of the same size/type) in one step - -ltid_abcrpsd: multivariable LTI sysID with - -25.02.09 - -modified mint_v2i to remove -repeated rows in its argument - - -20.02.09 - -modified msspoly/mtimes.m -to improve efficiency and memory handling - -new: msspoly/nnewton.m -Newton iterations with normalization - -new: msspoly/mono0.m -A new monomial generation function - -new: msspoly/mono_down.m -A new monomial generation function - -new: nlid_siso* function set -Nonlinear SISO system identification \ No newline at end of file diff --git a/spot_chk.m b/spot_chk.m deleted file mode 100644 index c0de285..0000000 --- a/spot_chk.m +++ /dev/null @@ -1,53 +0,0 @@ -% check a spot distribution - -fprintf(' Testing SPOT ...\n') - -fprintf('\n Test with mss_test1...') -h=mss_test1; -if abs(h+0.002519287852965)>1e-5, - error(' FAILED'); -else - fprintf(' OK\n') -end - -fprintf('\n Test with mss_test2...') -h=mss_test2; -if abs(h-2)>1e-5, - error(' FAILED'); -else - fprintf(' OK\n') -end - -fprintf('\n Test with mss_test3...') -h=mss_test3; -if abs(h-0.013647876833795)>1e-5, - error(' FAILED'); -else - fprintf(' OK\n') -end - -fprintf('\n Test with ltid_uy2ab_test...') -h=ltid_uy2ab_test; -if abs(h(2)-0.290068077153242)>1e-5, - error(' FAILED'); -else - fprintf(' OK\n') -end - -fprintf('\n Test with ltid_vw2ab_psd_test...') -h=ltid_vw2ab_psd_test(3,2,100,-1); -if h>0.031, - error(' FAILED'); -else - fprintf(' OK\n') -end - -fprintf('\n Test with nlid_fl_test1...') -h=nlid_fl_test1; -if h>0.05, - error(' FAILED'); -else - fprintf(' OK\n') -end - -fprintf('\n Test OK\n') \ No newline at end of file diff --git a/spot_install.m b/spot_install.m index 7bc078e..59e9644 100644 --- a/spot_install.m +++ b/spot_install.m @@ -12,9 +12,11 @@ addpath([potdir s 'util']); addpath([potdir s 'mint']); addpath([potdir s 'mss']); +addpath([potdir s 'spotopt']); +addpath([potdir s 'spotopt/solvers']); fprintf('\n compiling the binaries...') cd('bin'); mex mss_gset.c mex mss_gsum.c cd('..'); -fprintf('\n Done.\n') \ No newline at end of file +fprintf('\n Done.\n') diff --git a/spot_manual.pdf b/spot_manual.pdf deleted file mode 100644 index 435df28..0000000 Binary files a/spot_manual.pdf and /dev/null differ diff --git a/spot_manual.tex b/spot_manual.tex deleted file mode 100644 index af73b51..0000000 --- a/spot_manual.tex +++ /dev/null @@ -1,388 +0,0 @@ - -\documentclass[12pt]{article} -\usepackage{amsmath} % assumes amsmath package installed -\usepackage{amssymb} % assumes amsmath package installed -\usepackage{amsfonts} -\usepackage{graphicx} - -\title{ -SPOT\\ (Systems Polynomial Optimization Tools)\\ Manual -} - - -\author{Alexandre Megretski$\ ^\dagger$% -\thanks{$\ ^\dagger$LIDS, EECS, MIT, - {\tt\small ameg@mit.edu}}% -} - - -\newtheorem{thm}{Theorem} -\newtheorem{prop}{Proposition} -\newtheorem{cor}{Corollary} -\newtheorem{lem}{Lemma} -\newcounter{remark} -\newcounter{example} - -\newenvironment{rmk}{\addtocounter{remark}{1}\smallbreak\noindent - {\bf Remark \theremark}}{\hfill$\Box$\smallbreak} -\newenvironment{exmp}{\refstepcounter{example}\smallbreak\noindent - {\bf Example \theexample .}}{\hfill$\Box$\smallbreak} -\newenvironment{dfn}{\smallbreak\noindent{\bf -Definition.} }{\hfill$\Box$\smallbreak} -\newenvironment{pf}{\smallbreak\noindent{\it Proof. }}{\hfill$\Box$\smallbreak} -\newenvironment{pf*}[1]{\smallbreak\noindent{\it #1}}{\hfill$\Box$\smallbreak} - -\makeatother -\ifx\Box\undefined - \makeatletter\input{oldlfont.sty} \makeatother -\fi - - -\newcommand{\bb}[1]{\mbox{{\boldmath$#1$\unboldmath}}} -\newcommand{\bbo}[1]{\mbox{\rm{\bf #1}}} -\newcommand{\cs}{\mbox{\rm const}\cdot} -\newcommand{\rf}[1]{(\ref{#1})} -\newcommand{\cl}[1]{{\cal #1}} -\newcommand{\hh}[1]{\breve{#1}} -\newcommand{\col}{\mbox{\rm col}} -\newcommand{\tx}[1]{\mbox{\rm{#1}}} -\newcommand{\tre}[1]{2\mbox{\rm{Re}}\{#1\}} -\newcommand{\rre}[1]{\mbox{\rm{Re}}\{#1\}} -\def\Box{\hbox{\hskip 1pt \vrule width 8pt height 6pt depth 1.5pt - \hskip 1pt}} - -\renewcommand{\a}{\alpha} -\renewcommand{\b}{\beta} -\newcommand{\e}{\epsilon} -\newcommand{\g}{\gamma} -\newcommand{\Ga}{\Gamma} -\renewcommand{\d}{\delta} -\newcommand{\D}{\Delta} -\newcommand{\ka}{\kappa} -\newcommand{\w}{\omega} -\newcommand{\G}{\Gamma} - -\newcommand{\Wo}{\Omega} -\newcommand{\W}{\Omega} -\renewcommand{\r}{\rho} -\renewcommand{\t}{\theta} -\newcommand{\Th}{\Theta} -\newcommand{\s}{\sigma} -\renewcommand{\S}{\Sigma} -\newcommand{\n}{\nabla} -\renewcommand{\L}{\Lambda} - -\newcommand{\EQ}[2]{\begin{equation}\label{#1}#2\end{equation}} -\newcommand{\AR}[2]{\left[\begin{array}{#1}#2\end{array}\right]} -\newcommand{\ARR}[2]{\begin{array}{#1}#2\end{array}} -\newcommand{\OR}[2]{\left\{\begin{array}{#1}#2\end{array}\right.} -\newcommand{\NI}{\noindent} -\newcommand{\SK}{\vskip5mm} -\newcommand{\HD}[2]{\SK\NI{\bf #1.}\ \ #2\SK} -\newcommand{\la}{\langle} -\newcommand{\ra}{\rangle} -\newcommand{\sat}{\mbox{{\rm sat}}} -\newcommand{\sgn}{\mbox{{\rm sgn}}} -\newcommand{\dzn}{\mbox{{\rm dzn}}} -\newcommand{\dom}{\mbox{{\rm Dom}}} -\newcommand{\wto}{\to^*} -\newcommand{\LL}[2]{L_{#1}^{#2}} -\renewcommand{\sp}[2]{\langle #1,#2\rangle} -\newcommand{\rank}{\mbox{{\rm rank{\ }}}} -\newcommand{\sys}[4]{\left(\begin{array}{c|c}#1\\ \hline #3 -\end{array}\right)} -\newcommand{\EPS}[4]{ -\begin{center} -\resizebox{#2cm}{#3cm}{\rotatebox{#4}{\includegraphics{#1.eps}}} -\end{center} -} - -\newcommand{\RR}{\mathbb{R}} -\newcommand{\ZZ}{\mathbb{Z}} - - - -\begin{document} - - - -\maketitle -\thispagestyle{empty} -\pagestyle{empty} - - -SPOT (Systems Polynomial Optimization Tools) -is a MATLAB toolbox written as an alternative -implementation of SOSTools to be used in -implementing a class of -nonlinear system identification algorithms. -It was tested with {\tt MATLAB 7.8.0} -(R2009a). -SPOT provides its own matrix multivariable polynomial variable -class {\tt msspoly} for handling elementary polynomial operations, -a special class {\tt mssprog} for defining convex optimization problems -(to be solved by {\tt SeDuMi}) in terms -of polynomial identities and self-dual cones, -and a set of -functions for identification of linear and nonlinear dynamical systems. - -\section{Installation} -POT is distributed in the form of -compressed archives {\tt spotDDMMYY.zip}, -where {\tt DDMMYY} indicates the date of release (for example, -{\tt pot110110.zip} was released on January 11, 2010). - -Create directory {\tt spot} and extract {\tt spotDDMMYY.zip} into it. -Start MATLAB, and run {\tt spot\_install.m} from the {\tt pot} -directory. The script sets the path for POT, and compiles some binaries: -\begin{verbatim} ->> spot_install - - Installing SPOT in C:\home\matlab\spot: - updating the path... - compiling the binaries... - Done. ->> -\end{verbatim} - -Once SPOT is installed, you can check that it works by running -{\tt spot\_chk.m}: -\begin{verbatim} ->> spot_chk -\end{verbatim} - - - - -\section{Multivariable Polynomials} -The {\tt @msspoly} environment handles matrix polynomials in multiple variables. -Individual variables in {\tt @msspoly} have identifiers which begin with a -{\sl single} character, which may be followed by a non-negative integer number. -While, in principle, a variable identifier can begin with -any MATLAB-recognized character, the only ones which are safe to use in applications -are the 52 Latin alphabet letters -\begin{verbatim} -ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz -\end{verbatim} -Other characters can be used for automated variable definitions. For example, -in the {\tt mssprog} environment (used to define convex optimization programs -in terms of polynomial equations which are linear with respect to the optimized -({\sl decision}) variables, but have arbitrary dependence on other ({\sl abstract}) -variables, -the identifiers beginning with a ``{\tt @}'' are reserved for hidden decision -variables, while the identifiers beginning with a ``{\tt \#}'' are reserved for -hidden abstract variables. - -\subsection{Defining {\tt @msspoly} Variables} -Use {\tt msspoly.m}: -\begin{verbatim} ->> f=msspoly('v') - -[ v ] -\end{verbatim} -defines {\tt xx} as an {\tt @msspoly} polynomial {\tt f}$=f: f(v)=v$. -{\tt g=msspoly('v',k)} where {\tt k} is a positive integer -will define a {\tt k}-by-1 vector of different variables, as in -\begin{verbatim} ->> g=msspoly('t',3) - -[ t0 ] -[ t1 ] -[ t2 ] -\end{verbatim} -Using {\tt g=msspoly('v',[a b])} prodices a vector of {\tt a} variables with indexing -starting with {\tt b}, as in -\begin{verbatim} ->> g=msspoly('A',[2 1]) - -[ A1 ] -[ A2 ] -\end{verbatim} - - -\subsection{``Free'' and ``Simple'' {\tt @msspoly} Variables} -An {\tt @msspoly} variable is called {\sl free} if it is a matrix of -independent scalar variables. An {\tt @msspoly} variable is called {\sl simple} -if it is a column of -independent scalar variables and constants. For example, in -\begin{verbatim} ->> f1=msspoly('x',7); ->> f2=f1(1:6); ->> f3=reshape(f1,3,2); ->> f4=[f2;f2;1]; ->> f5=f1*f1'; -\end{verbatim} -the resulting free variables are {\tt f1}, {\tt f2}, {\tt f3}, -the simple variables are -{\tt f1}, {\tt f2}, {\tt f3}, {\tt f4}, while {\tt f5} is not free and not simple. - - - -\subsection{Handling {\tt @msspoly} Variables} -A number of functions of the {\tt @msspoly} environment have the -standard meaning: -\begin{itemize} -\item {\tt ctranspose.m} (as in {\tt z=x'}) -\item {\tt horzcat.m} (as in {\tt z=[x y]}) -\item {\tt minus.m} (as in {\tt z=x-y}) -\item {\tt mtimes.m} (as in {\tt z=x*y}) -\item {\tt plus.m} (as in {\tt z=x+y}) -\item {\tt uminus.m} (as in {\tt z=-x}) -\item {\tt uplus.m} (as in {\tt z=+x}) -\item {\tt vertcat.m} (as in {\tt z=[x;y]}) -\item {\tt subsasgn.m} (as in {\tt x(2)=y}) -\item {\tt reshape.m} -\item {\tt isempty.m} -\item {\tt isscalar.m} -\item {\tt length.m} -\item {\tt repmat.m} -\item {\tt size.m} -\item {\tt sum.m} -\end{itemize} -Other functions are close to their expected definitions, with minor -modifications or restrictions: -\begin{itemize} -\item {\tt decomp.m}: decomposes an {\tt @msspoly} variable into a vector of its -free variables, and matrices of degrees and coefficients of its terms; -\item {\tt deg.m}: gives the a single number degree; can be used with a second -argument, which must be a {\sl free} {\tt @msspoly} variable, in which case the -degree with respect to the independent variables listed in the second argument -is computed; -\item {\tt diag.m}: produces a diagonal {\tt @msspoly} matrix when the input is a row -or a column; otherwise extracts the diagonal as a column vector; -\item {\tt diff.m}: the second (required) and third (optional) arguments -must be free {\tt @msspoly} -variables; with two arguments, the first argument must be a column, and -the result is the matrix of the partial derivatives of the first -argument with respect to the second; with three -arguments, the first argument can have arbitrary dimensions, -and the result is the derivative of the first -argument with respect to the second in the direction provided by the third; -\item {\tt double.m}: converts a constant {\tt @msspoly} to {\tt double}, -otherwise returns character '?'; -\item {\tt isfunction.m}: the second argument must be a free {\tt @msspoly}; -true iff the firts argument is a function of the second; -\item {\tt mono.m}: produces column vector of all monomials from the argument; -\item {\tt mpower.m} (as in {\tt z=x\^{}y}): argument -{\tt y} must be a non-negative integer; -\item {\tt mrdivide.m} (as in {\tt z=x/y}): argument -{\tt y} must be a non-singular {\tt double}; -\item {\tt newton.m}: applies Newton method iterations to try to solve -approximately systems of polynomial equations; -\item {\tt recomp.m}: the inverse of {\tt decomp.m}; -\item {\tt subs.m}: a restricted substitution routine, allows to replace, -in the first argument, -the independent variables from the second argument -(must be a {\sl free} {\tt @msspoly}) -by the corresponding components of the third argument (must be a -{\sl simple} {\tt @msspoly}); -\item {\tt subsref.m} (as in {\tt z=x(1:2)} or {\tt z=x.n}): for the first type -of call, works as expected; for the second, {\tt x.m} and {\tt x.n} return -the dimensions, while {\tt x.s} returns the internal {\tt @msspoly} structure -of {\tt x} (something that only a developer of new {\tt @msspoly} code would need); -\item {\tt trace.m}: the usual sum of the diagonal elements, -but non-square arguments are admissible, too. -\end{itemize} - -\section{MSS Programs} - -MSS stands for ``Modified Sums of Squares''\footnote{or for ``Meager Sums of Squares'', -``Magnificent Sums of Squares'', etc., just not ``Alternative Sums of Squares'', -though that's what it is. } -The {\tt @mssprog} environment allows its user to define -matrix decision variables ranging over certain convex sets which are, in {\tt SeDuMi} -terminology, self-dual cones, to impose linear constraints in terms of polynomial -identities, to call {\tt SeDuMi} to optimize the decision variables, and to -extract the resulting optimal values. - -\subsection{{\tt @mssprog} Operations} -To initialize a blanc MSS program, use {\tt mssprog.m}, as in -\begin{verbatim} -pr=mssprog; -\end{verbatim} -The most straightforward way of adding items to an MSS program is by using the -{\tt '.'} subsassignments: -\begin{itemize} -\item {\tt pr.free=x} registers the elements of {\tt x} as free ({\tt SeDuMi}-style) decision -variables ({\tt x} must be a {\sl free} {\tt @msspoly}); -\item {\tt pr.pos=x} registers the elements of {\tt x} as positive ({\tt SeDuMi}-style) decision -variables ({\tt x} must be a {\sl free} {\tt @msspoly}); -\item {\tt pr.lor=x} registers the columns of {\tt x} as Lorentz cone ({\tt SeDuMi}-style) decision -variables ({\tt x} must be a {\sl free} {\tt @msspoly} with at least two rows); -\item {\tt pr.rlor=x} registers the columns of {\tt x} as rotated Lorentz cone -({\tt SeDuMi}-style) decision -variables ({\tt x} must be a {\sl free} {\tt @msspoly} with at least three rows); -\item {\tt pr.psd=x} registers the elements of every column of -{\tt x} as the components of positive semidefinite ({\tt SeDuMi}-style) decision -variables ({\tt x} must be a {\sl free} {\tt @msspoly} with {\tt nchoosek(m+1,2)} -rows to generate an {\tt m}-by-{\tt m} symmetric matrix: use -{\tt y=mss\_v2s(x(:,k)} to re-shape the {\tt k}-th column -of {\tt x} into the corresponding symmetric matrix); -\item {\tt pr.eq=x} registers equality {\tt x==0} with MSS program {\tt pr} -({\tt x} must be an {\tt @msspoly} which is linear with respect to the vector -of all independent variables which are registered with {\tt pr} as decision -parameters); -\item {\tt pr.sos=x} registers the constraint that -all scalar components of x must be sums of squares of polynomials -({\tt x} must be an {\tt @msspoly} which is linear with respect to the vector -of all independent variables which are registered with {\tt pr} as decision -parameters); -\item {\tt pr.sss=x} registers the constraint that -all {\tt u'*x*u}, where \newline -{\tt u=msspoly('\#',size(x,1))}, must be a sum of squares -({\tt x} must be a square-sized {\tt @msspoly} which is linear -with respect to the vector -of all independent variables which are registered with {\tt pr} as decision -parameters); -\item {\tt pr.sedumi=r} calls {\tt SeDuMi} to find the values of the decision -variables which minimize {\tt r} ({\tt r} must be a scalar {\tt @msspoly} -which is a linear function of the decision parameters). -\end{itemize} -To extract the optimized polynomials, use the {\tt '()'} subsreferencing: -\begin{itemize} -\item {\tt y=pr(x)}: {\tt y} is the result of substituting the optimized values -of decision variables into {\tt @msspoly} x; -\item {\tt y=pr(\{x\})}: same as {\tt double(pr(x))}. -\end{itemize} - -For example, the following code (contained in {\tt mss\_test3.m}) -finds the minimal value of $r$ for which -the polynomial $4x^4y^6+rx^2-xy^2+y^2$ is a sum of squares: -\begin{verbatim} -x=msspoly('x'); % define the variables -y=msspoly('y'); -r=msspoly('r'); -q=4*(x^4)*(y^6)+r*(x^2)-x*(y^2)+y^2; % the SOS polynomial -pr=mssprog; % initialize MSS program -pr.free=r; % register r as free -pr.sos=q; % register sos constraint -pr.sedumi=r; % minimize r -pr({r}) % get the optimal r -\end{verbatim} - - -\section{System Identification} -Functions from the {\tt nlid} directory are designed to help solving -nonlinear system identification problems. -Functions from the {\tt nlid} directory treat the linear time invariant (LTI) -case. - -\subsection{Identification of Symmetric Passive Transfer Matrices} -Function {\tt ltid\_passive.m} is a user interface to a number of -algorithms to aid in converting frequency samples of a -marginally stable symmetric transfer -matrix (such as impedance of a passive circuit) -to a reduced order state space model. - - - -\subsection{Nonlinear System Identification} -Currently, the following examples from the {\tt nlid} directory appear to work: -\begin{itemize} -\item {\tt nlid\_io\_lti\_old\_test1.m}: LTI system identification with -POT; -\item {\tt nlid\_miso0\_test1.m}: memoryless NL system identification; -\item {\tt nlid\_fl\_test1.m}: more powerful memoryless NL system identification. -\end{itemize} - -\end{document} diff --git a/@SPOTSOSProg/SPOTSOSProg.m b/spotopt/@SPOTSOSProg.old/SPOTSOSProg.m similarity index 100% rename from @SPOTSOSProg/SPOTSOSProg.m rename to spotopt/@SPOTSOSProg.old/SPOTSOSProg.m diff --git a/@SPOTSOSProg/test_lb.m b/spotopt/@SPOTSOSProg.old/test_lb.m similarity index 100% rename from @SPOTSOSProg/test_lb.m rename to spotopt/@SPOTSOSProg.old/test_lb.m diff --git a/@SPOTSOSSoln/SPOTSOSSoln.m b/spotopt/@SPOTSOSSoln/SPOTSOSSoln.m similarity index 100% rename from @SPOTSOSSoln/SPOTSOSSoln.m rename to spotopt/@SPOTSOSSoln/SPOTSOSSoln.m diff --git a/@SPOTSQLProg/SPOTSQLProg.m b/spotopt/@SPOTSQLProg/SPOTSQLProg.m similarity index 100% rename from @SPOTSQLProg/SPOTSQLProg.m rename to spotopt/@SPOTSQLProg/SPOTSQLProg.m diff --git a/@SPOTSQLProg/test_nl_fitting.m b/spotopt/@SPOTSQLProg/test_nl_fitting.m similarity index 100% rename from @SPOTSQLProg/test_nl_fitting.m rename to spotopt/@SPOTSQLProg/test_nl_fitting.m diff --git a/@SPOTSQLProg/test_nl_fitting_dual.m b/spotopt/@SPOTSQLProg/test_nl_fitting_dual.m similarity index 100% rename from @SPOTSQLProg/test_nl_fitting_dual.m rename to spotopt/@SPOTSQLProg/test_nl_fitting_dual.m diff --git a/@SPOTSQLProg/test_nl_fitting_mssprog.m b/spotopt/@SPOTSQLProg/test_nl_fitting_mssprog.m similarity index 100% rename from @SPOTSQLProg/test_nl_fitting_mssprog.m rename to spotopt/@SPOTSQLProg/test_nl_fitting_mssprog.m diff --git a/@SPOTSQLProg/test_nl_fitting_primal.m b/spotopt/@SPOTSQLProg/test_nl_fitting_primal.m similarity index 100% rename from @SPOTSQLProg/test_nl_fitting_primal.m rename to spotopt/@SPOTSQLProg/test_nl_fitting_primal.m diff --git a/@SPOTSQLProg/test_nl_fitting_primal_foo.m b/spotopt/@SPOTSQLProg/test_nl_fitting_primal_foo.m similarity index 100% rename from @SPOTSQLProg/test_nl_fitting_primal_foo.m rename to spotopt/@SPOTSQLProg/test_nl_fitting_primal_foo.m diff --git a/@SPOTSQLProg/test_nl_fitting_setup.m b/spotopt/@SPOTSQLProg/test_nl_fitting_setup.m similarity index 100% rename from @SPOTSQLProg/test_nl_fitting_setup.m rename to spotopt/@SPOTSQLProg/test_nl_fitting_setup.m diff --git a/@SPOTSQLSoln/SPOTSQLSoln.m b/spotopt/@SPOTSQLSoln/SPOTSQLSoln.m similarity index 100% rename from @SPOTSQLSoln/SPOTSQLSoln.m rename to spotopt/@SPOTSQLSoln/SPOTSQLSoln.m diff --git a/spotopt/@SpotProgPreProc/SpotProgPreProc.m b/spotopt/@SpotProgPreProc/SpotProgPreProc.m new file mode 100644 index 0000000..bf95a61 --- /dev/null +++ b/spotopt/@SpotProgPreProc/SpotProgPreProc.m @@ -0,0 +1,5 @@ +classdef SpotProgPreProc + methods (Abstract) + [sdpout,G,h,log] = preProcess(sdpin); + end +end \ No newline at end of file diff --git a/spotopt/@spotprog/spotprog.m b/spotopt/@spotprog/spotprog.m new file mode 100644 index 0000000..480e034 --- /dev/null +++ b/spotopt/@spotprog/spotprog.m @@ -0,0 +1,157 @@ +classdef spotprog + properties + longname = ''; + sdp=[]; % Current SDP Problem Description + + numVar = 0; % Number of variables. + + % Mapping taking sdp variables into user facing variables. + variableExpr = []; + preProc = {}; + log = {}; + end + + methods ( Access = private ) + function nm = varName(prog) + nm = [prog.sdp.name 'var']; + end + + function [prog,v] = newVariables(prog,var_sdp) + prog.variableExpr = [ prog.variableExpr ; var_sdp ]; + v = msspoly(prog.varName,[length(var_sdp) prog.numVar]); + prog.numVar = prog.numVar + length(var_sdp); + end + + function sdpVar = toSDPVariables(pr,var) + sdpVar = subs(var,pr.variables,pr.variableExpr); + end + + end + + methods + function prog = spotprog(shortname,longname) + if nargin >= 1, + prog.sdp = spotsdp(shortname); + else + prog.sdp = spotsdp; + end + + if nargin >= 2, + prog.longname = longname; + else + prog.longname = 'default'; + end + end + + function v = variables(pr) + v = msspoly(pr.varName,pr.numVar); + end + + function pr = withEqs(pr,e) + pr.sdp = pr.sdp.withEqs(pr.toSDPVariables(e)); + end + function pr = withPos(pr,e) + pr.sdp = pr.sdp.withPos(pr.toSDPVariables(e)); + end + function pr = withLor(pr,e) + pr.sdp = pr.sdp.withLor(pr.toSDPVariables(e)); + end + function pr = withRLor(pr,e) + pr.sdp = pr.sdp.withRLor(pr.toSDPVariables(e)); + end + function pr = withPSD(pr,e) + pr.sdp = pr.sdp.withPSD(pr.toSDPVariables(e)); + end + function pr = withBlkPSD(pr,e) + pr.sdp = pr.sdp.withBlkPSD(pr.toSDPVariables(e)); + end + + + function [pr,v] = newFree(pr,dim) + [pr.sdp,s] = pr.sdp.newFree(dim); + [pr,v] = newVariables(pr,s); + end + + function [pr,v] = newPos(pr,dim) + [pr.sdp,s] = pr.sdp.newPos(dim); + [pr,v] = newVariables(pr,s); + end + + function [pr,v] = newLor(pr,dim) + [pr.sdp,s] = pr.sdp.newLor(dim); + [pr,v] = newVariables(pr,s(:)); + v = reshape(v,size(s)); + end + + function [pr,v] = newRLor(pr,dim) + [pr.sdp,s] = pr.sdp.newRLor(dim); + [pr,v] = newVariables(pr,s(:)); + v = reshape(v,size(s)); + end + + function [pr,v] = newPSD(pr,dim) + [pr.sdp,s] = pr.sdp.newPSD(dim); + [pr,v] = newVariables(pr,mss_s2v(s)); + v = mss_v2s(v); + end + + + function [pr,v] = newBlkPSD(pr,dim) + [pr.sdp,s] = pr.sdp.newBlkPSD(dim); + [pr,v] = newVariables(pr,s(:)); + v = reshape(v,size(s)); + end + + function pr = setPreProc(pr,preproc) + % + % pr = addPreProc(pr,preproc) + % + % pr -- a spotprog + % preproc -- a SpotProgPreProc object. + % + pr.preProc{end+1} = preproc; + end + + + function [sdp,A,b,obj] = prepare(pr,objective) + objective = subs(objective,pr.variables,pr.variableExpr); + + [A,b] = spot_decomp_linear(pr.variableExpr,pr.sdp.variables); + b = -b; + + % A*(Ai*z+bi)+b + sdp = pr.sdp; + for i = 1:length(pr.preProc) + [sdp,Ai,bi,log{i}] = pr.preProc{i}.preProcess(sdp); + A = A*Ai; + b = b+A*bi; + end + + obj = subs(objective,pr.variables,b+A*sdp.variables); + end + + function sol = minimize(pr,solver,objective) + % + % sol = pr.minimize(objective,solver) + % [sol,err] = pr.minimize(objective,solver) + % + % objective -- 1-by-1 msspoly linear in pr.variables + % solver -- spotsolver + % + % sol -- If the given solver can handle the problem type, sol is + % a spotsol. If the program cannot be handled + % err = 1 for two output arguments and an error is + % thrown o.w. + ; + % First solve the SDP + [sdp,A,b,obj] = prepare(pr,objective); + + int_sol = solver.minimize(sdp,obj); + + sol = spotsdpsol(pr,int_sol.info,pr.variables,... + int_sol.primalSolution,A,b); + end + + end + +end diff --git a/spotopt/@spotsdp/TODO b/spotopt/@spotsdp/TODO new file mode 100644 index 0000000..deed703 --- /dev/null +++ b/spotopt/@spotsdp/TODO @@ -0,0 +1,2 @@ +add withRLor +switch dual representation. \ No newline at end of file diff --git a/spotopt/@spotsdp/approxSDPbySDD.m b/spotopt/@spotsdp/approxSDPbySDD.m new file mode 100644 index 0000000..9f47fde --- /dev/null +++ b/spotopt/@spotsdp/approxSDPbySDD.m @@ -0,0 +1,56 @@ +function [approx,oldVariables,newVariables] = approxSDPbySDD(pr) +% +% [approx,oldVar,newExpr] = approxSDPbySDD(pr) +% +% pr -- spotsdp. +% oldVar -- Variables replaced. +% newExpr -- Expression of old variables in new variables. +% approx -- New spotsdp. +% +% Approx is a new SDP with each sdp constraint and variable +% replaced by a scaled diagonally dominant constraint. + + % Step one: Construct new program with new cone variables + approx = spotsdp; + + fold = pr.freeVariables; + [approx,fnew] = approx.newFree(pr.numFree); + + pold = pr.posVariables; + [approx,pnew] = approx.newPos(pr.numPos); + + lold = pr.lorVariables; + [approx,lnew] = approx.newLors(pr.lorDim); + + rold = pr.rlorVariables; + [approx,rnew] = approx.newRLors(pr.rlorDim); + + psdold = pr.psdVariables; + psddim = pr.psdDim; + psdnew = []; + for i = 1:length(pr.psdDim) + [approx,psdApprox] = newSDD(approx,pr.psdDim(i)); + psdnew = [ psdnew ; psdApprox ]; + end + + oldVariables = [ fold ; pold ; lold ; rold ; psdold ]; + newVariables = [ fnew ; pnew ; lnew ; rnew ; psdnew ]; + + % Step Two: Walk through constraints, replacing new variables. + + [pcstr,dim] = pr.posConstraints(); + approx = approx.withPos(subs(pcstr,oldVariables,newVariables)); + + [lcstr,dim] = pr.lorConstraints(); + approx = approx.withLor(subs(lcstr,oldVariables,newVariables),dim); + + [rcstr,dim] = pr.rlorConstraints(); + approx = approx.withRLor(subs(rcstr,oldVariables,newVariables),dim); + + [~,dim] = pr.psdConstraints(); + for i = 1:length(dim) + scstr = subs(pr.psdConstraint(i),oldVariables,newVariables); + [approx,slack] = newSDD(approx,pr.psdCstrDim(i)); + approx = approx.withEqs(scstr - slack); + end +end \ No newline at end of file diff --git a/spotopt/@spotsdp/newDD.m b/spotopt/@spotsdp/newDD.m new file mode 100644 index 0000000..325928d --- /dev/null +++ b/spotopt/@spotsdp/newDD.m @@ -0,0 +1,43 @@ +function [prog,psdApprox] = newDD(prog,n) +% +% [prout,psd] = pr.newSDD(n) +% +% Adds a "scaled diagonally dominant matrix variable. +% O.w. known as a factor-width 2 PSD matrix. +% + M = nchoosek(n,2); + [prog,d] = prog.newPos(n); + [prog,p] = prog.newPos(M); + [prog,m] = prog.newPos(M); + + psdP = [ 1 1 1 ]'*p'; + psdM = [ 1 -1 1 ]'*m'; + + + N2 = mss_s2v(reshape(1:n^2,n,n)); + + % Enumerate unique pairs (i,j), i < j. + [i,j]=ind2sub([n n],N2); + ondiag = i == j; + id = i(ondiag); + jd = i(ondiag); + + offdiag = i ~= j; + i = i(offdiag); + j = j(offdiag); + + % Produce row/column that each 2-by-2 should be mapped to. + I = [ i' ; i' ; j' ]; + J = [ i' ; j' ; j' ]; + + % Produce the matrix index that corresponds to these. + Ind = sub2ind([n n],I(:),J(:)); + VInd = mss_match(N2,Ind); + S = sparse(VInd,(1:length(VInd))',ones(length(VInd),1),length(N2),length(VInd)); + + Indd = sub2ind([n n],id,jd); + VIndd = mss_match(N2,Indd); + Sd = sparse(VIndd,(1:length(VIndd))',ones(length(VIndd),1),length(N2),length(VIndd)); + + psdApprox = S*psdP(:) + S*psdM(:) + Sd*d; +end \ No newline at end of file diff --git a/spotopt/@spotsdp/newSDD.m b/spotopt/@spotsdp/newSDD.m new file mode 100644 index 0000000..6bf4842 --- /dev/null +++ b/spotopt/@spotsdp/newSDD.m @@ -0,0 +1,33 @@ +function [prog,psdApprox] = newSDD(prog,n) +% +% [prout,psd] = pr.newSDD(n) +% +% Adds a "scaled diagonally dominant matrix variable. +% O.w. known as a factor-width 2 PSD matrix. +% + M = nchoosek(n,2); + [prog,r] = prog.newRLor([3 M]); + + % Map to 2-by-2 SDP matrices in upper triangular format. + psd2 = [ sqrt(2) 0 0 ; 0 0 1 ; 0 sqrt(2) 0]*r; + + N2 = mss_s2v(reshape(1:n^2,n,n)); + + % Enumerate unique pairs (i,j), i < j. + [i,j]=ind2sub([n n],N2); + offdiag = i ~= j; + i = i(offdiag); + j = j(offdiag); + + % Produce row/column that each 2-by-2 should be mapped to. + I = [ i' ; i' ; j' ]; + J = [ i' ; j' ; j' ]; + + % Produce the matrix index that corresponds to these. + Ind = sub2ind([n n],I(:),J(:)); + % Now map matrix indices into the upper triangular vector format. + VInd = mss_match(N2,Ind); + S = sparse(VInd,(1:length(VInd))',ones(length(VInd),1),length(N2),length(VInd)); + + psdApprox = S*psd2(:); +end \ No newline at end of file diff --git a/spotopt/@spotsdp/primalize.m b/spotopt/@spotsdp/primalize.m new file mode 100644 index 0000000..4dadf0f --- /dev/null +++ b/spotopt/@spotsdp/primalize.m @@ -0,0 +1,50 @@ +function [prout,G,h] = primalize(prin) +% +% [prout,G,h] = primalize(prin) +% +% Converts a program into the standard primal with free +% variables form via the introduction of slack variables. +% +% The matrices G,h are constructed so that +% +% prg.variables = G*spPrg.variables + h. +% + prout = prin; + prout.posCnst = []; + prout.lorCstr = []; + prout.lorCstrDim = []; + prout.rlorCstr = []; + prout.rlorCstrDim = []; + prout.psdCstr = []; + prout.psdCstrDim = []; + + [pos,dim] = prin.posConstraints(); + if ~isempty(dim) && dim ~= 0 + [prout,slack] = prout.newPos(dim); + prout = prout.withEqs(pos-slack); + end + + [lor,dim] = prin.lorConstraints(); + if ~isempty(dim) + [prout,slack] = prout.newLors(dim); + prout = prout.withEqs(lor-slack); + end + + [rlor,dim] = prin.rlorConstraints(); + if ~isempty(dim) + [prout,slack] = prout.newRLors(dim); + prout = prout.withEqs(rlor-slack); + end + + [psd,dim] = prin.psdConstraints(); + if ~isempty(dim) + [prout,slack] = prout.newPSDs(dim); + prout = prout.withEqs(psd-slack); + end + + h = zeros(size(prin.variables)); + + [var,pow,Coeff] = decomp(prout.variables); + mtch = match(var,prin.variables); + G = Coeff(:,mtch)'; +end \ No newline at end of file diff --git a/spotopt/@spotsdp/spotsdp.m b/spotopt/@spotsdp/spotsdp.m new file mode 100644 index 0000000..e481586 --- /dev/null +++ b/spotopt/@spotsdp/spotsdp.m @@ -0,0 +1,547 @@ +classdef spotsdp + properties + % Developer's Notes about Internal Representation: + % + % The program consists of a collection of variable dimensions: + % + % name -- Character prefix for variables from the program. + % + % + % posNum -- 1-by-1 Positive integer (number of non-negative variables). + % freeNum -- 1-by-1 Positive integer (number of non-negative variables). + % psdDim -- Npsd-by-1 array of positive integers. Each represents + % psdDim(i)-by-psdDim(i) dim. variable. + % lorDim -- Nlor-by-1 array of positive integers. Size of + % Lorentz cones (n indicates x(1)^2 >= + % sum_i=2^n x(i)^2 ) + % rlorDim -- Nrlor-by-1 array of positive integers, similar + % for rotated Lorentz cones. + % + % Variables are named @psdi, @lori, @posi @rlri, where @ is + % replaced by 'name' and 'i' is a running counter. + % + % + % + name = '@'; + + posNum = 0; + freeNum = 0; + psdDim = []; + lorDim = []; + rlorDim = []; + + psdCstr = []; + psdCstrDim = []; + + lorCstr = []; + lorCstrDim = []; + rlorCstr = []; + rlorCstrDim = []; + posCnst = []; + + equations = []; + end + + methods (Static) + function n = psdDimToNo(d) + n=(d+1).*d/2; + end + + function [d,v] = psdNoToDim(n) + % + % 0 = d^2 + d - 2n + % (sqrt(1 + 8n) - 1)/2 + % + d=round((sqrt(1+8*n)-1)/2); + if spotsdp.psdDimToNo(d) ~= n + d = NaN; + v = 0; + else + v = 1; + end + end + + end + + + methods ( Access = private ) +% These private functions define the names of variables +% generated by the program. + + function nm = freeName(pr) + nm = [pr.name 'fr']; + end + function nm = posName(pr) + nm = [pr.name 'pos']; + end + function nm = psdName(pr) + nm = [pr.name 'psd']; + end + function nm = lorName(pr) + nm = [pr.name 'lor']; + end + function nm = rlorName(pr) + nm = [pr.name 'rlr']; + end + + function flag = legalEq(pr,eq) + if ~isa(eq,'msspoly') + flag = 0; + else + flag = realLinearInDec(pr,eq); + end + end + end + + % ---- General modeling methods (i.e. posing inequalities) + methods + function pr=spotsdp(name) + % pr=spotsdp(prefix) + % + % prefix -- Scalar character, legal name for msspoly. + % + % Returns: + % pr -- New program, decision variables begin with + % the character prefix. + % + % + % spotsdp objects model SDP/SOCP/LP cone programs. + % + % The feasible set of these programs is represented in a + % mixture of standard primal and standard dual form: + % + % (F) x in K1, y free, + % A1.x + A2.y = b, + % D1.x + D2.y + e in K2, + % + % where K1 and K2 are products of the SDP, SOCP and LP + % cones of various dimensions. + % + % After the feasible set has been constructed, an + % optimization problem can be solved via prg.optimze(): + % + % minimize c'x + f'y subj. to. (F) + % + % Solution of these problems by primal dual solvers such + % as SeDuMi or SDPT3 requires converting the problem to + % standard primal or dual form. The choice of which + % conversion to apply can be forced as an optimizatio parameter. + % + % + % + if nargin > 0 + if ~ischar(name) || length(name) > 1 + error('Program name must be a scalar character.'); + else + msspoly(name); + end + pr.name = name; + end + end + + function flag = realLinearInDec(pr,exp) + [x,pow,Coeff] = decomp(exp); + [~,xid] = isfree(x); + [~,vid] = isfree(pr.variables); + flag = ~(any(mss_match(vid,xid) == 0) | ... + any(pow(:) > 1) | ... + any(imag(Coeff(:)) ~= 0)); + end + + % Generate variables of a given type. + function f = freeVariables(pr) + f = msspoly(pr.freeName,pr.numFree); + end + function p = posVariables(pr) + p = msspoly(pr.posName,pr.numPos); + end + function [l,dim] = lorVariables(pr) + l = msspoly(pr.lorName,pr.numLor); + dim = pr.lorDim; + end + + function [r,dim] = rlorVariables(pr) + r = msspoly(pr.rlorName,pr.numRLor); + dim = pr.rlorDim; + end + + function [p,dim] = psdVariables(pr) + p = msspoly(pr.psdName,pr.numPSD); + dim = pr.psdDim; + end + + function p = psdVariable(pr,i) + if ~spot_isIntGE(i,1) || i > length(pr.psdDim) + error(['PSD variable ' num2str(i) ' does not exist.']); + end + i0 = sum(spotsdp.psdDimToNo(pr.psdDim(1:i-1))); + n = spotsdp.psdDimToNo(pr.psdDim(i)); + p = mss_v2s(msspoly(pr.psdName,[n,i0])); + end + + function [pcstr,dim] = posConstraints(pr) + pcstr = pr.posCnst; + dim = length(pcstr); + end + + function [lcstr,dim] = lorConstraints(pr) + lcstr = pr.lorCstr; + dim = pr.lorCstrDim; + end + + function [rcstr,dim] = rlorConstraints(pr) + rcstr = pr.rlorCstr; + dim = pr.rlorCstrDim; + end + + + function [psdcstr,dim] = psdConstraints(pr) + dim = pr.psdCstrDim; + psdcstr = pr.psdCstr; + end + + function [scstr,dim] = psdConstraint(pr,i) + if ~spot_hasSize(i,[1 1]) || ~spot_isIntGE(i,1) || ... + i > length(pr.psdCstrDim) + error('Bad Index.'); + end + + i0 = sum(spotsdp.psdDimToNo(pr.psdCstrDim(1:i-1))); + dim = pr.psdCstrDim(i); + scstr = pr.psdCstr(i0+(1:spotsdp.psdDimToNo(dim))); + end + + + function v = variables(pr) + % v = variables(pr) + % v -- msspoly column of primal variables for the program pr. + v = [ pr.freeVariables + pr.posVariables + pr.lorVariables + pr.rlorVariables + pr.psdVariables]; + end + + function n = numPos(pr) + n = pr.posNum; + end + function n = numFree(pr) + n = pr.freeNum; + end + function n = numPSD(pr) + n = sum(spotsdp.psdDimToNo(pr.psdDim)); + end + function n = numLor(pr) + n = sum(pr.lorDim); + end + function n = numRLor(pr) + n = sum(pr.rlorDim); + end + + function m = numEq(pr) + m = length(pr.equations); + end + + + function [pr,Q] = newPSD(pr,dim) + if ~spot_hasSize(dim,[1 1]) || ~spot_isIntGE(dim,1) + error('Dimension must be scalar positive integer.'); + end + n = spotsdp.psdDimToNo(dim); + + Q = mss_v2s(msspoly(pr.psdName,[n pr.numPSD])); + + pr.psdDim = [pr.psdDim dim]; + end + + function [pr,Q] = newPSDs(pr,dim) + if isempty(dim), Q = []; return; end + + if size(dim,1) ~= 1 || ~spot_isIntGE(dim,1) + error('Dimension must be scalar positive integers.'); + end + + n = sum(spotsdp.psdDimToNo(dim)); + + Q = msspoly(pr.psdName,[n pr.numPSD]); + + pr.psdDim = [ pr.psdDim dim ]; + end + + function [pr,Qs] = newBlkPSD(pr,dim) + if ~spot_hasSize(dim,[1 2]) || ~spot_isIntGE(dim,1) + error('Dimension must be 1x2 positive integer.'); + end + + n = spotsdp.psdDimToNo(dim(1)); + + Qs = reshape(msspoly(pr.psdName,[n*dim(2) pr.numPSD]),n,dim(2)); + pr.psdDim = [pr.psdDim dim(1)*ones(1,dim(2))]; + end + + function [pr,p] = newPos(pr,dim) + if ~spot_hasSize(dim,[1 1]) || ~spot_isIntGE(dim,0) + error('Dimension must be scalar positive integer.'); + end + + if dim == 0, p = []; + else + p = msspoly(pr.posName,[dim pr.numPos]); + pr.posNum = pr.posNum+dim; + end + end + + function [pr,f] = newFree(pr,dim) + if spot_hasSize(dim,[1 1]) + dim = [ dim 1]; + end + if ~spot_hasSize(dim,[1 2]) || ~spot_isIntGE(dim,1) + error('Dimension must be 1-by-1 or 1-by-2 positive integer.'); + end + + f = reshape(msspoly(pr.freeName,[prod(dim) pr.numFree]),dim); + + pr.freeNum = pr.freeNum+prod(dim); + end + + function [pr,l] = newLor(pr,dim) + if spot_hasSize(dim,[1 1]), dim = [dim 1]; end + if ~spot_hasSize(dim,[1 2]) || ~spot_isIntGE(dim,1) + error('Dimension must be 1x2 positive integer.'); + end + + l = reshape(msspoly(pr.lorName,[prod(dim) pr.numLor]),dim(1),dim(2)); + + pr.lorDim = [pr.lorDim dim(1)*ones(1,dim(2))]; + end + + function [pr,l] = newLors(pr,dim) + if isempty(dim), l = []; return; end + + if size(dim,1) ~= 1 || ~spot_isIntGE(dim,1) + error('Dimension must be 1xn positive integer.'); + end + + l = msspoly(pr.lorName,[sum(dim) pr.numLor]); + pr.lorDim = [ pr.lorDim dim ]; + end + + function [pr,r] = newRLor(pr,dim) + if spot_hasSize(dim,[1 1]), dim = [dim 1]; end + if ~spot_hasSize(dim,[1 2]) || ~spot_isIntGE(dim,1) + error('Dimension must be 1x2 positive integer.'); + end + + r = reshape(msspoly(pr.rlorName,[prod(dim) pr.numRLor]),dim(1),dim(2)); + + pr.rlorDim = [pr.rlorDim dim(1)*ones(1,dim(2))]; + end + + function [pr,r] = newRLors(pr,dim) + if isempty(dim), r = []; return; end + + if size(dim,1) ~= 1 || ~spot_isIntGE(dim,2) + error('Dimension must be 1xn integer, >= 2.'); + end + + r = msspoly(pr.rlorName,[sum(dim) pr.numRLor]); + pr.rlorDim = [ pr.rlorDim dim ]; + end + + function [pr] = withEqs(pr,eq) + if ~pr.legalEq(eq) + error(['Equations must be an msspoly linear in ' ... + 'decision parameters.']); + end + eq = eq(:); + + pr.equations = [pr.equations ; eq]; + end + + %-- + function [pr] = withPos(pr,exp) + if ~isa(exp,'msspoly') + error('Argument must be a column of msspoly expressions.'); + end + exp = exp(:); + + pr.posCnst = [pr.posCnst ; exp]; + end + + function [pr] = withPSD(pr,exp,dim) + if ~isa(exp,'msspoly') + error('Argument must be an msspoly.'); + end + + % It is non-square. + if size(exp,1) ~= size(exp,2) || nargin > 2, + if size(exp,2) > 1, + error('Argument must be square or a column.'); + end + if nargin < 2, + error(['For column argument, dim must be specified']); + end + + if sum(spotsdp.psdDimToNo(dim)) ~= length(exp) + error(['Dimension does not agree with length of ' ... + 'column argument.']); + end + + pr.psdCstr = [pr.psdCstr ; exp ]; + pr.psdCstrDim = [ pr.psdCstrDim dim]; + else % Non-square and dim not specified. + dim = size(exp,1); + exp = mss_s2v(exp); + pr = pr.withPSD(exp,dim); + end + end + + function [pr] = withBlkPSD(pr,exp) + if ~isa(exp,'msspoly') + error('Argument must be an msspoly.'); + end + [dim,v] = spotsdp.psdNoToDim(size(exp,1)); + if ~v + error('Argument wrong size.'); + end + + pr = pr.withPSD(exp(:),dim*ones(1,size(exp,2))); + end + + function [pr] = withLor(pr,exp,dim) + if nargin > 2 && isempty(dim), + return; + end + if ~isa(exp,'msspoly') + error('Argument must be an msspoly.'); + end + + if nargin == 2, + pr = pr.withLor(exp(:),size(exp,1)*ones(1,size(exp,2))); + else + if size(dim,1) ~= 1 || ~spot_isIntGE(dim,1), + error('Dimensions must be row of positive integers.'); + end + + if sum(dim) ~= length(exp), + error(['Dimensions and expression length do not match.']); + end + pr.lorCstr = [pr.lorCstr ; exp]; + pr.lorCstrDim = [ pr.lorCstrDim dim ]; + end + end + + function [pr] = withRLor(pr,exp,dim) + if nargin > 2 && isempty(dim), + return; + end + if ~isa(exp,'msspoly') + error('Argument must be an msspoly.'); + end + + if nargin == 2, + pr = pr.withLor(exp(:),size(exp,1)*ones(1,size(exp,2))); + else + if size(dim,1) ~= 1 || ~spot_isIntGE(dim,2), + error('Dimensions must be row of integers, >= 2.'); + end + + if sum(dim) ~= length(exp), + error(['Dimensions and expression length do not match.']); + end + + pr.rlorCstr = [pr.rlorCstr ; exp]; + pr.rlorCstrDim = [ pr.rlorCstrDim dim ]; + end + end + + end + + methods + function f = isStandardDualForm(prg) + f = isempty(prg.equations) && ... + prg.numPos == 0 && ... + prg.numPSD == 0 && ... + prg.numLor == 0 && ... + prg.numRLor == 0; + end + +% Transform programs into standard forms (maybe move this out?) + + function [prgout,G,h] = standardDual(prg) + % + % [spPrg,G,h] = standardDual(prg) + % + % Converts a program into the standard dual form. + % + % Conic variables are replaced by new free variables. + % + % Then, a lower dimensional parameterization: + % + % [ x; y] = Gz + h + % + % is found with [A1 A2] h = b, [A1 A2] G = 0. + % + % + + function prgout = moveConstraints(prgout,prg,free) + if ~isempty(prg.posCnst) + prgout = prgout.withPos(subs(prg.posCnst,prg.variables,free)); + end + end + + + function prgout = removeConic(prg) + prgout = spotsdp(prg.name); + [prgout,free] = prgout.newFree(length(prg.variables)); + + if prg.numPos > 0 + mtch = match(prg.variables,prg.posVariables); + prgout = prgout.withPos(free(mtch)); + end + + if prg.numLor > 0 + error('Did not support Lorentz cone yet.'); + end + if prg.numRLor > 0 + error('Did not support Rotated Lorentz cone yet.'); + end + + prgout = moveConstraints(prgout,prg,free); + prgout = prgout.withEqs(subs(prg.equations,prg.variables,free)); + end + + function prgout = removeEquality(prg) + % Now resolve equality constraints. + if length(prg.equations) == 0, + prgout = prg; + G = speye(prg.numFree); + h = sparse([],[],[],prg.numFree,1); + return; + end + + prgout = spotsdp(prg.name); + + [A,b] = spot_decomp_linear(prg.equations,prg.variables); + + % TODO: Return to this setting to preserve sparsity. + + h = A\b; + + [Q,R] = qr(A'); + n = max(find(sum(abs(R),2))); + G = Q(:,n+1:end); % New basis + + [prgout,z] = prgout.newFree(size(G,2)); + + prgout = moveConstraints(prgout,prg,G*z+h); + end + + prgout = removeConic(prg); + prgout = removeEquality(prgout); + + end + end + + +end \ No newline at end of file diff --git a/spotopt/@spotsdp/toDual.m b/spotopt/@spotsdp/toDual.m new file mode 100644 index 0000000..7e2cf9d --- /dev/null +++ b/spotopt/@spotsdp/toDual.m @@ -0,0 +1,120 @@ +function [dl,dobj] = toDual(pr,pobj) +% +% [dl,dobj] = toDual(pr,pobj) +% +% pr -- spotsdp +% pobj -- 1-by-1 linear msspoly, function of pr.variables. +% dual -- spotsdp +% dobj -- 1-by-1 msspoly +% +% +% pr is given by the mixed primal dual description: +% +% f free, x in K1, A1*x + A2*f = b, F1*x + F2*f - g in K2 +% +% pobj = + (to be minimized) +% +% dl is the dual conic program derived from the Lagrangian: +% +% + - + - +% +% i.e.: +% +% y free, c1'+A1'*y-F1'*s in K1*, s in K2*, A2'*y-F2'*s + c2' =0 +% +% with dobj = - . + ; + f = pr.freeVariables; + x = [ pr.posVariables ; pr.lorVariables ; pr.rlorVariables ; pr.psdVariables]; + + % x in K1, determine inner product for this cone. + Q1 = inner_product(length(pr.posVariables)+... + length(pr.lorVariables)+... + length(pr.rlorVariables),pr.psdDim); + + + dl = spotsdp; + + if isempty(pr.equations) + A1Ty = zeros(length(x),1); + A2Ty = zeros(length(f),1); + bTy = 0; + else + [A,b] = spot_decomp_linear(pr.equations,[x;f]); + A1 = A(:,1:length(x)); A2 = A(:,length(x)+1:end); + + % Construct dual with variables. + [dl,y] = dl.newFree(length(b)); + + A2Ty = A2'*y; + A1Ty = Q1\A1'*y; + bTy = b'*y; + end + + [c,~] = spot_decomp_linear(pobj,[x;f]); + c1 = c(:,1:length(x)); c2 = c(:,length(x)+1:end); + c1T = c1'; + + % Next, determine the cone K2* + [pcstr,pdim] = pr.posConstraints(); + [lcstr,ldim] = pr.lorConstraints(); + [rcstr,rdim] = pr.rlorConstraints(); + [scstr,sdim] = pr.psdConstraints(); + + [dl,pdl] = dl.newPos(pdim); + [dl,ldl] = dl.newLors(ldim); + [dl,rdl] = dl.newRLors(rdim); + [dl,sdl] = dl.newPSDs(sdim); + + s = [ pdl ; ldl ; rdl ; sdl ]; + + % Determine Equations + K2 = [ pcstr ; lcstr ; rcstr ; scstr]; + + % Need to determine the inner-product for this cone. + % Particularly the PSD cone. + Q2 = inner_product(length(pcstr)+length(lcstr)+length(rcstr),sdim); + + if ~isempty(K2) + [F,g] = spot_decomp_linear(K2,[x;f]); + F1 = F(:,1:length(x)); F2 = F(:,length(x)+1:end); + + F1Ts = Q1\F1'*Q2*s; + F2Ts = F2'*Q2*s; + gTs = g'*Q2*s; + else + F1Ts = zeros(size(c1T)); + F2Ts = zeros(size(c2')); + gTs = 0; + end + + dl = dl.withEqs(A2Ty-F2Ts+c2'); + + % Finally determine cone K1* + K1star = c1T+A1Ty-F1Ts; + + ipos = (1:pr.numPos)'; + ilor = pr.numPos+(1:pr.numLor)'; + irlor = pr.numPos+pr.numLor+(1:pr.numRLor)'; + ipsd = pr.numPos+pr.numLor+pr.numRLor+(1:pr.numPSD)'; + + dl = dl.withPos(K1star(ipos)); + dl = dl.withLor(K1star(ilor),pr.lorDim); + dl = dl.withRLor(K1star(irlor),pr.rlorDim); + dl = dl.withPSD(K1star(ipsd),pr.psdDim); + + dobj = gTs-bTy; + + function Q = inner_product(n,dim) + I = 2*ones(sum(spotsdp.psdDimToNo(dim)),1); + + i0 = 0; + for i = 1:length(dim) + ii = cumsum(1:dim(i)); + I(i0+ii) = I(i0+ii) - 1; + i0 = i0 + spotsdp.psdDimToNo(dim(i)); + end + + Q = blkdiag(speye(n),spdiags(I,0,length(I),length(I))); + end +end \ No newline at end of file diff --git a/@spotsqlsol/SPOTSQLSoln.m b/spotopt/@spotsdpsol/SPOTSQLSoln.m similarity index 100% rename from @spotsqlsol/SPOTSQLSoln.m rename to spotopt/@spotsdpsol/SPOTSQLSoln.m diff --git a/spotopt/@spotsdpsol/spotsdpsol.m b/spotopt/@spotsdpsol/spotsdpsol.m new file mode 100644 index 0000000..899395b --- /dev/null +++ b/spotopt/@spotsdpsol/spotsdpsol.m @@ -0,0 +1,41 @@ +classdef spotsdpsol + properties + program = []; + user_variables = []; + primalSolution = []; + G = []; + h = []; + info = []; + end + + methods + function f = primalInfeasible(sol) + f = sol.info.pinf; + end + function f = dualInfeasible(sol) + f = sol.info.dinf; + end + + function sol = spotsdpsol(prog,info,uservars,psol,G,h) + if nargin < 5, G = speye(length(psol)); end + if nargin < 6, h = sparse(length(psol),1); end + + sol.program = prog; + sol.info = info; + sol.user_variables = uservars; + sol.primalSolution = psol; + sol.G = G; + sol.h = h; + + end + + function ev = eval(sol,exp) + if sol.primalInfeasible + error('Primal Infeasible, cannot evalutate primal variables.'); + end + + ev = subs(exp,[sol.user_variables],... + [sol.G*sol.primalSolution+sol.h]); + end + end +end \ No newline at end of file diff --git a/@spotsosprg/spotsosprg.m b/spotopt/@spotsosprg/spotsosprg.m similarity index 100% rename from @spotsosprg/spotsosprg.m rename to spotopt/@spotsosprg/spotsosprg.m diff --git a/spotopt/@spotsosprog/spotsosprog.m b/spotopt/@spotsosprog/spotsosprog.m new file mode 100644 index 0000000..1035a4e --- /dev/null +++ b/spotopt/@spotsosprog/spotsosprog.m @@ -0,0 +1,200 @@ +classdef spotsosprog < spotprog + properties + sosExpr = []; + end + methods (Access = private) + function [flag,indet] = realPolyLinearInDec(pr,exp) + [x,pow,Coeff] = decomp(exp); + + mtch = match(x,pr.variables); + + flag = ~(any(any(pow(:,mtch(mtch~=0)) > 1)) | ... + any(imag(Coeff(:)) ~= 0)); + + if ~flag, + indet = []; + else + mtch = match(pr.variables,x); + indet = x(find(mtch==0)); + if any(istrig(indet)) + flag = 0; + indet = []; + end + end + end + + function [flag,tIn,pIn] = trigPolyLinearInDec(pr,expr) + [x,pow,Coeff] = decomp(expr); + + mtch = match(x,pr.variables); + + % TODO: conj. symmetric check. + flag = ~(any(any(pow(:,mtch(mtch~=0)) > 1))); + + if ~flag, + indet = []; + else + mtch = match(pr.variables,x); + indet = x(find(mtch==0)); + msk = istrig(indet); + tIn = indet(find(msk)); + pIn = indet(find(~msk)); + end + end + end + + methods (Access = private, Static) + function phi = buildGramBasis(expr,decvar) + if ~spot_hasSize(expr,[1 1]) + error('buildGramBasis expects a scalar polynomial.'); + end + + [var,pow,M] = decomp(expr); + mtch = match(var,decvar); + b = 1:length(var); + b(mtch(mtch ~= 0)) = []; + indet = var(b); + + if length(indet) == 0 + phi = 1; + return; + end + + pow = pow(:,b); + + exponent_m = spot_build_gram_basis(pow); + + phi = recomp(indet,exponent_m,eye(size(exponent_m,1))); + end + end + + methods (Access = protected) + function [pr,Q,phi,y,basis] = buildSOSDecompPrimal(pr,expr) + if ~spot_hasSize(expr,[1 1]) + error('buildSOSDecomp expects a scalar polynomial.'); + end + + decvar = pr.variables; + + phi = spotsosprog.buildGramBasis(expr,decvar); + + [pr,Q] = pr.newPSD(length(phi)); + + decvar = [decvar ; mss_s2v(Q)]; + sosCnst = expr-phi'*Q*phi; + + A = diff(sosCnst,decvar); + b = subs(sosCnst,decvar,0*decvar); + [var,pow,Coeff] = decomp([b A].'); + + %[pr,y] = pr.withEqs(Coeff'*[1;decvar]); + [pr] = pr.withEqs(Coeff'*[1;decvar]); + y = []; + basis = recomp(var,pow,eye(size(pow,1))); + end + end + + methods + function pr = spotsosprog(varargin) + pr@spotprog(varargin{:}); + end + + + function [pr,poly,coeff] = newSOSPoly(pr,basis,n) + if nargin < 3, n = 1; end + [pr,poly,coeff] = newFreePoly(pr,basis,n); + pr = pr.withSOS(poly); + end + + function [pr] = withSOS(pr,expr) + if ~pr.realPolyLinearInDec(expr) + error(['Coefficients must be real, indeterminates ' ... + 'non-trigonometric, and expression must ' ... + 'be linear in decision variables.']); + end + tokens = length(pr.sosExpr) + (1:prod(size(expr))); + pr.sosExpr = [ pr.sosExpr ; expr(:)]; + end + + function [pr] = withPolyEqs(pr,expr) + if ~pr.realPolyLinearInDec(expr) + error(['Coefficients must be real, indeterminates ' ... + 'non-trigonometric, and expression must ' ... + 'be linear in decision variables.']); + end + + expr = expr(:); + decvar = pr.variables; + + [indet,pow,M] = decomp(expr,decvar); + + monom = recomp(indet,pow,eye(size(pow,1))); + + [I,J,S] = find(M); + + [pr,y] = pr.withEqs(S); + + basis = monom(J); + end + + function [pr] = withSOSMatrix(pr,expr) + [lindec,indet] = pr.realPolyLinearInDec(expr); + if ~lindec + error(['Coefficients must be real, indeterminates ' ... + 'non-trigonometric, and expression must ' ... + 'be linear in decision variables.']); + end + + if size(expr,1) ~= size(expr,2) || size(expr,1) == 0 + error('Expression must be a square non-empty matrix.'); + end + + x = msspoly('x',size(expr,1)); + expr = anonymize(expr,'y',indet); + [pr,tokens] = withSOS(pr,x'*expr*x); + end + + function [pr] = withUniTrigSOSMatrix(pr,expr) + if size(expr,1) ~= size(expr,2) || size(expr,1) == 0 + error('Expression must be a square non-empty matrix.'); + end + + [lindec,trigIndet,polyIndet] = pr.trigPolyLinearInDec(expr); + if ~lindec + error(['Expression ' ... + 'must be linear in decision variables.']); + elseif ~isempty(polyIndet) || (length(trigIndet) ~= 1) + error('Only a single, trigonometric indeterminate allowed.'); + end + + + % How to handle these? Special case? Maybe... + + end + + function n = numSOS(pr) + n = length(pr.sosExpr); + end + + function [pr,poly,coeff] = newFreePoly(pr,basis,n) + if nargin < 3, n = 1; end + [pr,coeff] = pr.newFree(length(basis)*n); + poly = reshape(coeff,n,length(basis))*basis; + end + + + function sol = minimize(pr,varargin) + if nargin < 2, objective = msspoly(0); end + + Q = cell(pr.numSOS,1); + phi = cell(pr.numSOS,1); + y = cell(pr.numSOS,1); + basis = cell(pr.numSOS,1); + for i = 1:pr.numSOS + [pr,Q{i},phi{i},y{i},basis{i}] = pr.buildSOSDecompPrimal(pr.sosExpr(i)); + end + + sol = minimize@spotprog(pr,varargin{:}); + end + end +end \ No newline at end of file diff --git a/@spotsossol/spotsossol.m b/spotopt/@spotsossol/spotsossol.m similarity index 100% rename from @spotsossol/spotsossol.m rename to spotopt/@spotsossol/spotsossol.m diff --git a/@spotsqlprg/SPOTSQLProg.m b/spotopt/@spotsqlprg/SPOTSQLProg.m similarity index 100% rename from @spotsqlprg/SPOTSQLProg.m rename to spotopt/@spotsqlprg/SPOTSQLProg.m diff --git a/@spotsqlprg/spotsqlprg.m b/spotopt/@spotsqlprg/spotsqlprg.m similarity index 100% rename from @spotsqlprg/spotsqlprg.m rename to spotopt/@spotsqlprg/spotsqlprg.m diff --git a/@spotsqlprg/test_fitting.m b/spotopt/@spotsqlprg/test_fitting.m similarity index 100% rename from @spotsqlprg/test_fitting.m rename to spotopt/@spotsqlprg/test_fitting.m diff --git a/@spotsqlprg/test_nl_fitting.m b/spotopt/@spotsqlprg/test_nl_fitting.m similarity index 100% rename from @spotsqlprg/test_nl_fitting.m rename to spotopt/@spotsqlprg/test_nl_fitting.m diff --git a/@spotsqlprg/test_rqp.m b/spotopt/@spotsqlprg/test_rqp.m similarity index 100% rename from @spotsqlprg/test_rqp.m rename to spotopt/@spotsqlprg/test_rqp.m diff --git a/spotopt/TODO b/spotopt/TODO new file mode 100644 index 0000000..b1918e7 --- /dev/null +++ b/spotopt/TODO @@ -0,0 +1,20 @@ +Overall TODO: +------------- + +Finish spotprog: + - Basic functionality: layer around spotsdp. + - Add "constraints" which optionally can have + primal form and dual form representations. + - Write SOSConstraint class. + +Improve spotsdpsolver: + - spotsdpsol shoudl store the solver as well. + - Maybe add an "advice" functionality like sysid toolbox? + - options interfaces + - Need a interface for deciding can it be solved. + - Add the ability to "preprocess" + + + + + diff --git a/spotopt/solvers/@spotsdpsolver/spotsdpsolver.m b/spotopt/solvers/@spotsdpsolver/spotsdpsolver.m new file mode 100644 index 0000000..ab9487a --- /dev/null +++ b/spotopt/solvers/@spotsdpsolver/spotsdpsolver.m @@ -0,0 +1,166 @@ +classdef spotsdpsolver + methods + function can = canSolve(solvr,sdp) + can = 0; + end + + function sol = minimizePrimalForm(solvr,sdp,objective) + error('minimizePrimalForm not implemented.'); + end + + function sol = minimizeDualForm(solvr,sdp,objective) + error('minimizeDualForm not implemented.'); + end + + function sol = minimize(solver,pr,objective) + + if pr.isStandardDualForm() + sol = solver.minimizeDualForm(pr,objective); + else + sol = solver.minimizePrimalForm(pr,objective); + end + end + end + + + + methods (Static) + function [As,bs] = linearToSedumi(lin,vall,varNo,KvarCnt) + [A,bs] = spot_decomp_linear(lin,vall); + [i,j,s] = find(A); + As = sparse(i,varNo(j),s,size(A,1),KvarCnt); + end + + function [psdVarNo,psdVarNoSymm] = upperTriToFullVarNo(npsd,psdDim) + % Assign column numbers to v. + psdVarNo = zeros(1,npsd); + psdVarNoSymm = zeros(1,npsd); + psdVarOff = 0; % Progress in variables, storing + % upper triangle. + psdRedVarOff = 0; % Progress in variables, storing + % entire matrix. + for i = 1:length(psdDim) + n = psdDim(i); + m = n*(n+1)/2; + psdVarNo(psdVarOff + (1:m)) = psdRedVarOff+mss_s2v(reshape(1:n^2,n,n)); + psdVarNoSymm(psdVarOff + (1:m)) = psdRedVarOff+mss_s2v(reshape(1:n^2,n,n)'); + psdVarOff = psdVarOff + m; + psdRedVarOff = psdRedVarOff + n^2; + end + end + + function [A,b,c,K,G,h,varNo] = sdpToPrimalSedumiFormat(pr,objective) + objective = msspoly(objective); + if ~realLinearInDec(pr,objective) + error('Objective must be real and linear in dec. variables.'); + end + + user_variables = pr.variables; + + [pr,G,h] = pr.primalize(); + + % First, construct structure with counts of SeDuMi + % variables. + K = struct(); + K.f = pr.freeNum; + K.l = pr.posNum; + K.q = pr.lorDim; + K.r = pr.rlorDim; + K.s = pr.psdDim; + + KvarCnt = K.f+K.l+sum(K.q)+sum(K.r)+sum(K.s.^2); + + + v = [ pr.freeVariables + pr.posVariables + pr.lorVariables + pr.rlorVariables ]; + + vpsd = pr.psdVariables; + + vall = [v;vpsd]; + + [psdVarNo] = spotsdpsolver.upperTriToFullVarNo(length(vpsd),pr.psdDim); + + varNo = [ 1:length(v) length(v)+psdVarNo]; + + [A,b] = spotsdpsolver.linearToSedumi(pr.equations,vall,varNo,KvarCnt); + [c,~] = spotsdpsolver.linearToSedumi(objective,vall,varNo,KvarCnt); + end + + function [A,b,c,K,G,h] = sdpToDualSedumiFormat(pr,objective) + user_variables = pr.variables; + + objective = msspoly(objective); + if ~realLinearInDec(pr,objective) + error('Objective must be real and linear in dec. variables.'); + end + + if ~pr.isStandardDualForm() + error(['Right now the program must be in standard dual ' ... + 'format.']); + + [pr,G,h] = pr.standardDual(); + else + G = speye(size(user_variables,1)); + h = sparse([],[],[],size(user_variables,1),size(user_variables,2)); + end + + objective = subs(objective,user_variables,G*pr.variables+h); + + % + % maximize b'y + % c-A'y in K. + % + % (1) Decide on the appropriate cone sizes. + % (2) Construct A, c. + + % First, construct structure with counts of SeDuMi + % variables. + K = struct(); + K.f = 0; + + [pos,K.l] = pr.posConstraints(); + [lor,K.q] = lorConstraints(pr); + [rlor,K.r] = rlorConstraints(pr); + + v = [ pos ; lor ; rlor ]; + + [vpsd,K.s] = pr.psdConstraints(); + + KvarCnt = K.f+K.l+sum(K.q)+sum(K.r)+sum(K.s.^2); + + [psdVarNo,psdVarNoSymm] = spotsdpsolver.upperTriToFullVarNo(length(vpsd),K.s); + + varNo = [ 1:length(v) length(v)+psdVarNo]; + varNoSymm = [ 1:length(v) length(v)+psdVarNoSymm]; + varNoDiag = varNo(find(varNo == varNoSymm)); + + + % I need to isolate the term: c - A'y in K where K is + % the cone of the above dimensions. + + vall = [ v ; vpsd ]; + + [mAT,cm] = spot_decomp_linear(vall,pr.variables); + [b,~] = spot_decomp_linear(objective,pr.variables); + b = -b.'; + + [i,j,s] = find(mAT); + + mAT = sparse(varNo(i),j,s,KvarCnt,size(mAT,2)); + keep = varNo(i) ~= varNoSymm(i); + mAT = mAT+sparse(varNoSymm(i(keep)),... + j(keep),s(keep),KvarCnt,size(mAT,2)); + + [i,j,s] = find(-cm); + c = sparse(varNo(i),j,s,KvarCnt,1); + keep = varNo(i) ~= varNoSymm(i); + c = c + sparse(varNoSymm(i(keep)),j(keep),s(keep),KvarCnt,1); + + A = -mAT.'; + end + + + end +end \ No newline at end of file diff --git a/spotopt/solvers/@spotsolversdpnal/spotsolversdpnal.m b/spotopt/solvers/@spotsolversdpnal/spotsolversdpnal.m new file mode 100644 index 0000000..aaf8f51 --- /dev/null +++ b/spotopt/solvers/@spotsolversdpnal/spotsolversdpnal.m @@ -0,0 +1,44 @@ +classdef spotsolversdpnal < spotsdpsolver + properties + options = struct(); + end + + methods + + function solvr = spotsolversdpnal(options) + ; + % TODO: check that all options are legal. + if nargin > 0 + solvr.options = options; + else + solvr.options = struct(); + end + end + + function can = canSolve(solvr,sdp) + can = 1; + end + + function sol = minimize(solver,pr,objective) + error('Optimization in primal form not supported.'); + end + + function sol = minimizeDualForm(solver,pr,objective) + [A,b,c,K,G,h] = spotsdpsolver.sdpToDualSedumiFormat(pr,objective); + + [blk,Asdpt,Csdpt,bsdpt] = read_sedumi_sdpnal_0(A,full(b),c,K); + + [obj,X,y,Z,infosdpt] = sdpnal(blk,Asdpt,Csdpt,bsdpt); + info = struct('pinf',infosdpt.termcode == 1,... + 'dinf',infosdpt.termcode == 2); + + if info.dinf || info.pinf, + primalSol = NaN*ones(size(pr.variables,1),1); + else + primalSol = y; + end + + sol = spotsdpsol(pr,info,pr.variables,G*primalSol+h); + end + end +end diff --git a/spotopt/solvers/@spotsolversdpt3_4/spotsolversdpt3_4.m b/spotopt/solvers/@spotsolversdpt3_4/spotsolversdpt3_4.m new file mode 100644 index 0000000..e6f1152 --- /dev/null +++ b/spotopt/solvers/@spotsolversdpt3_4/spotsolversdpt3_4.m @@ -0,0 +1,40 @@ +classdef spotsolversdpt3_4 < spotsdpsolver + properties + options = struct(); + end + + methods + + function solvr = spotsolversdpt3_4(options) + ; + % TODO: check that all options are legal. + if nargin > 0 + solvr.options = options; + else + solvr.options = struct(); + end + end + + function can = canSolve(solvr,sdp) + can = sdp.isStandardDualForm(); + end + + function sol = minimizeDualForm(solver,pr,objective) + [A,b,c,K,G,h] = spotsdpsolver.sdpToDualSedumiFormat(pr,objective); + + [blk,Asdpt,Csdpt,bsdpt] = read_sedumi_sdpt3_4(A,full(b),c,K); + + [obj,X,y,Z,infosdpt] = sdpt3(blk,Asdpt,Csdpt,bsdpt); + info = struct('pinf',infosdpt.termcode == 1,... + 'dinf',infosdpt.termcode == 2); + + if info.dinf || info.pinf, + primalSol = NaN*ones(size(pr.variables,1),1); + else + primalSol = y; + end + + sol = spotsdpsol(pr,info,pr.variables,G*primalSol+h); + end + end +end diff --git a/spotopt/solvers/@spotsolversedumi/spotsolversedumi.m b/spotopt/solvers/@spotsolversedumi/spotsolversedumi.m new file mode 100644 index 0000000..f450cba --- /dev/null +++ b/spotopt/solvers/@spotsolversedumi/spotsolversedumi.m @@ -0,0 +1,60 @@ +classdef spotsolversedumi < spotsdpsolver + properties + options = struct(); + end + + + methods + + function solvr = spotsolversedumi(options) + ; + % TODO: check that all options are legal. + if nargin > 0 + solvr.options = options; + else + solvr.options = struct('fid',0); + end + end + + + function can = canSolve(solvr,sdp) + can = 1; + end + + + + function sol = minimizePrimalForm(solver,pr,objective) + if nargin < 2, objective = 0; end + + [A,b,c,K,G,h,varNo] = spotsdpsolver.sdpToPrimalSedumiFormat(pr,objective); + + [x,y,info] = sedumi(A,b,c,K,solver.options); + + if info.pinf, + primalSol = NaN*ones(size(varNo,2),1); + else + primalSol = x(varNo); + end + + primalSol = G*primalSol + h; + + sol = spotsdpsol(pr,info,pr.variables,primalSol); + end + + function sol = minimizeDualForm(solver,pr,objective) + + [A,b,c,K,G,h] = spotsdpsolver.sdpToDualSedumiFormat(pr,objective); + + [x,y,info] = sedumi(A,b,c,K,solver.options); + + + if info.dinf || info.pinf, + primalSol = NaN*ones(size(pr.variables,1),1); + else + primalSol = y; + end + + sol = spotsdpsol(pr,info,pr.variables,G*primalSol+h); + end + end +end diff --git a/spotopt/solvers/read_sedumi_sdpnal_0.m b/spotopt/solvers/read_sedumi_sdpnal_0.m new file mode 100644 index 0000000..9f5bdae --- /dev/null +++ b/spotopt/solvers/read_sedumi_sdpnal_0.m @@ -0,0 +1,188 @@ +%%******************************************************************* +%% Read in a problem in SeDuMi format. +%% +%% [blk,A,C,b,perm] = read_sedumi(fname,b,c,K) +%% +%% Input: fname.mat = name of the file containing SDP data in +%% SeDuMi format. +%% +%% Important note: Sedumi's notation for free variables "K.f" +%% is coded in SDPT3 as blk{p,1} = 'u', where +%% "u" is used for unrestricted variables. +%% +%% SDPNAL: +%% Copyright (c) 2008 by +%% Xinyuan Zhao, Defeng Sun, and Kim-Chuan Toh +%%****************************************************************** + + function [blk,Avec,C,b,perm] = read_sedumi(fname,b,c,K,smallblkdim); + + if (nargin < 5) + smallblkdim = 30; + end + + A = 0; + At = 0; + if isstr(fname) + %% + %% load the matlab file containing At, c, b, and K + %% + K.f = []; K.l = []; K.q = []; + compressed = 0; + if exist([fname,'.mat.gz']); + compressed = 1; + unix(['gunzip ', fname,'.mat.gz']); + elseif exist([fname,'.gz']); + compressed = 2; + unix(['gunzip ', fname,'.gz']); + elseif exist([fname,'.mat.Z']); + compressed = 3; + unix(['uncompress ', fname,'.mat.Z']); + elseif exist([fname,'.Z']); + compressed = 4; + unix(['uncompress ', fname,'.Z']); + end + if exist([fname,'.mat']) | exist(fname) + eval(['load ', fname]); + else + fprintf('*** Problem not found, please specify the correct path or problem. \n'); + blk = []; Avec = []; C = []; b = []; + return; + end + if (compressed == 1) + unix(['gzip ', fname,'.mat']); + elseif (compressed == 2) + unix(['gzip ', fname]); + elseif (compressed == 3) + unix(['compress ', fname,'.mat']); + elseif (compressed == 4) + unix(['compress ', fname]); + end + elseif (nargin < 4) + error('read_sedumi: need 4 input '); + else + A = fname; + end +%% + if exist('c','var') + if (size(c,1) == 1), c = c'; end; + end + if exist('C','var') + c = C; + if (size(c,1) == 1), c = c'; end; + end + if (size(b,1) == 1), b = b'; end; + if (norm(A,'fro') > 0) & (size(A,2) == length(b)); At = A; end +%% + if (norm(At,'fro')==0), At = A'; end; + [nn,mm] = size(At); if (max(size(c)) == 1); c = c*ones(nn,1); end; + if ~isfield(K,'f'); K.f = 0; end + if ~isfield(K,'l'); K.l = 0; end + if ~isfield(K,'q'); K.q = 0; end + if ~isfield(K,'s'); K.s = 0; end + if (K.f == 0) | isempty(K.f); K.f = 0; end; + if (K.l == 0) | isempty(K.l); K.l = 0; end; + if (sum(K.q) == 0) | isempty(K.q); K.q = 0; end + if (sum(K.s) == 0) | isempty(K.s); K.s = 0; end +%% +%% +%% + m = length(b); + rowidx = 0; idxblk = 0; + if ~(K.f == 0) + len = K.f; + idxblk = idxblk + 1; + blk{idxblk,1} = 'u'; blk{idxblk,2} = K.f; + Atmp = At(rowidx+[1:len],:); + Avec{idxblk,1} = Atmp; + C{idxblk,1} = c(rowidx+[1:len]); + perm{idxblk} = []; + rowidx = rowidx + len; + end + if ~(K.l == 0) + len = K.l; + idxblk = idxblk + 1; + blk{idxblk,1} = 'l'; blk{idxblk,2} = K.l; + Atmp = At(rowidx+[1:len],:); + Avec{idxblk,1} = Atmp; + C{idxblk,1} = c(rowidx+[1:len]); + perm{idxblk} = []; + rowidx = rowidx + len; + end + if ~(K.q == 0) + len = sum(K.q); + idxblk = idxblk + 1; + blk{idxblk,1} = 'q'; + if size(K.q,1) <= size(K.q,2); + blk{idxblk,2} = K.q; + else + blk{idxblk,2} = K.q'; + end + Atmp = At(rowidx+[1:len],:); + Avec{idxblk,1} = Atmp; + C{idxblk,1} = c(rowidx+[1:len]); + perm{idxblk} = []; + rowidx = rowidx + len; + end + if ~(K.s == 0) + blksize = K.s; + if (size(blksize,2) == 1); blksize = blksize'; end + blknnz = [0, cumsum(blksize.*blksize)]; + deblkidx = find(blksize > smallblkdim); + if ~isempty(deblkidx) + for p = 1:length(deblkidx) + idxblk = idxblk + 1; + n = blksize(deblkidx(p)); + pblk{1,1} = 's'; pblk{1,2} = n; + blk(idxblk,:) = pblk; + Atmp = At(rowidx+blknnz(deblkidx(p))+[1:n*n],:); + h = [1:n]'; e = ones(n,1); + tmp = triu(h*e' + e*((h-1).*h/2)'); + tmp = tmp(:); + tmp2 = sqrt(2)*triu(ones(n),1) + speye(n,n); + tmp2 = tmp2(:); + symidx = find(tmp(:)); + dd = tmp2(find(tmp2)); + n2 = n*(n+1)/2; + Avec{idxblk,1} = spdiags(dd,0,n2,n2)*Atmp(symidx,:); + Ctmp = c(rowidx+blknnz(deblkidx(p))+[1:n*n]); + Ctmp = mexmat(pblk,Ctmp,1); + C{idxblk,1} = 0.5*(Ctmp+Ctmp'); + perm{idxblk,1} = deblkidx(p); + end + end + spblkidx = find(blksize <= smallblkdim); + if ~isempty(spblkidx) + cnt = 0; cnt2 = 0; + spblksize = blksize(spblkidx); + nn = sum(spblksize.*spblksize); + nn2 = sum(spblksize.*(spblksize+1)/2); + pos = zeros(nn,1); + dd = zeros(nn2,1); + symidx = zeros(nn2,1); + for p = 1:length(spblkidx) + n = blksize(spblkidx(p)); + n2 = n*(n+1)/2; + pos(cnt+[1:n*n]) = rowidx+blknnz(spblkidx(p))+[1:n*n]; + h = [1:n]'; e = ones(n,1); + tmp = triu(h*e' + e*((h-1).*h/2)'); + tmp = tmp(:); + tmp2 = sqrt(2)*triu(ones(n),1) + speye(n,n); + tmp2 = tmp2(:); + symidx(cnt2+[1:n2]) = cnt+find(tmp(:)); + dd(cnt2+[1:n2]) = tmp2(find(tmp2)); + cnt = cnt + n*n; + cnt2 = cnt2 + n2; + end + idxblk = idxblk + 1; + blk{idxblk,1} = 's'; blk{idxblk,2} = blksize(spblkidx); + Atmp = At(pos,:); + Avec{idxblk,1} = spdiags(dd,0,length(dd),length(dd))*Atmp(symidx,:); + Ctmp = c(pos); + Ctmp = mexmat(blk(idxblk,:),Ctmp,1); + C{idxblk,1} = 0.5*(Ctmp+Ctmp'); + perm{idxblk,1} = spblkidx; + end + end +%% +%%******************************************************************* diff --git a/spotopt/solvers/read_sedumi_sdpt3_4.m b/spotopt/solvers/read_sedumi_sdpt3_4.m new file mode 100644 index 0000000..70107d6 --- /dev/null +++ b/spotopt/solvers/read_sedumi_sdpt3_4.m @@ -0,0 +1,222 @@ +%%******************************************************************* +%% Read in a problem in SeDuMi format. +%% +%% [blk,A,C,b,perm] = read_sedumi(fname,b,c,K) +%% +%% Input: fname.mat = name of the file containing SDP data in +%% SeDuMi format. +%% +%% Important note: Sedumi's notation for free variables "K.f" +%% is coded in SDPT3 as blk{p,1} = 'u', where +%% "u" is used for unrestricted variables. +%% +%% SDPT3: version 3.1 +%% Copyright (c) 1997 by +%% K.C. Toh, M.J. Todd, R.H. Tutuncu +%% Last Modified: 16 Sep 2004 +%%****************************************************************** + + function [blk,Avec,C,b,perm] = read_sedumi(fname,b,c,K,smallblkdim); + + if (nargin < 5) + smallblkdim = 50; + end + + A = 0; + At = 0; + if isstr(fname) + %% + %% load the matlab file containing At, c, b, and K + %% + K.f = []; K.l = []; K.q = []; + compressed = 0; + if exist([fname,'.mat.gz']); + compressed = 1; + unix(['gunzip ', fname,'.mat.gz']); + elseif exist([fname,'.gz']); + compressed = 2; + unix(['gunzip ', fname,'.gz']); + elseif exist([fname,'.mat.Z']); + compressed = 3; + unix(['uncompress ', fname,'.mat.Z']); + elseif exist([fname,'.Z']); + compressed = 4; + unix(['uncompress ', fname,'.Z']); + end + if exist([fname,'.mat']) | exist(fname) + eval(['load ', fname]); + else + fprintf('*** Problem not found, please specify the correct path or problem. \n'); + blk = []; Avec = []; C = []; b = []; + return; + end + if (compressed == 1) + unix(['gzip ', fname,'.mat']); + elseif (compressed == 2) + unix(['gzip ', fname]); + elseif (compressed == 3) + unix(['compress ', fname,'.mat']); + elseif (compressed == 4) + unix(['compress ', fname]); + end + elseif (nargin < 4) + error('read_sedumi: need 4 input '); + else + A = fname; + end +%% + if exist('c','var') + if (size(c,1) == 1), c = c'; end; + end + if exist('C','var') + c = C; + if (size(c,1) == 1), c = c'; end; + end + if (size(b,1) == 1), b = b'; end; + if (norm(A,'fro') > 0) & (size(A,2) == length(b)); At = A; end +%% + if (norm(At,'fro')==0), At = A'; end; + [nn,mm] = size(At); if (max(size(c)) == 1); c = c*ones(nn,1); end; + if ~isfield(K,'f'); K.f = 0; end + if ~isfield(K,'l'); K.l = 0; end + if ~isfield(K,'q'); K.q = 0; end + if ~isfield(K,'s'); K.s = 0; end + if (K.f == 0) | isempty(K.f); K.f = 0; end; + if (K.l == 0) | isempty(K.l); K.l = 0; end; + if (sum(K.q) == 0) | isempty(K.q); K.q = 0; end + if (sum(K.s) == 0) | isempty(K.s); K.s = 0; end +%% +%% +%% + m = length(b); + rowidx = 0; idxblk = 0; + if ~(K.f == 0) + len = K.f; + idxblk = idxblk + 1; + blk{idxblk,1} = 'u'; blk{idxblk,2} = K.f; + Atmp = At(rowidx+[1:len],:); + Avec{idxblk,1} = Atmp; + C{idxblk,1} = c(rowidx+[1:len]); + perm{idxblk} = []; + rowidx = rowidx + len; + end + if ~(K.l == 0) + len = K.l; + idxblk = idxblk + 1; + blk{idxblk,1} = 'l'; blk{idxblk,2} = K.l; + Atmp = At(rowidx+[1:len],:); + Avec{idxblk,1} = Atmp; + C{idxblk,1} = c(rowidx+[1:len]); + perm{idxblk} = []; + rowidx = rowidx + len; + end + if ~(K.q == 0) + len = sum(K.q); + idxblk = idxblk + 1; + blk{idxblk,1} = 'q'; + if size(K.q,1) <= size(K.q,2); + blk{idxblk,2} = K.q; + else + blk{idxblk,2} = K.q'; + end + Atmp = At(rowidx+[1:len],:); + Avec{idxblk,1} = Atmp; + C{idxblk,1} = c(rowidx+[1:len]); + perm{idxblk} = []; + rowidx = rowidx + len; + end + if ~(K.s == 0) + blksize = K.s; + if (size(blksize,2) == 1); blksize = blksize'; end + blknnz = [0, cumsum(blksize.*blksize)]; + deblkidx = find(blksize > smallblkdim); + if ~isempty(deblkidx) + for p = 1:length(deblkidx) + idxblk = idxblk + 1; + n = blksize(deblkidx(p)); + pblk{1,1} = 's'; pblk{1,2} = n; + blk(idxblk,:) = pblk; + Atmp = At(rowidx+blknnz(deblkidx(p))+[1:n*n],:); + %% + %% column-wise positions of upper triangular part + %% + tmp = triu(ones(n)); tmp = tmp(:); + idxtriu = find(tmp); + %% + %% row-wise positions of lower triangular part + %% + tmp = tril(reshape([1:n*n],n,n)); tmp = tmp(:); + idxtmp = find(tmp); + [dummy,idxsub] = sort(rem(tmp(idxtmp),n)); + idxtril = [idxtmp(idxsub(n+1:end));idxtmp(idxsub(1:n))]; + %% + tmp2 = sqrt(2)*triu(ones(n),1) + speye(n,n); + tmp2 = tmp2(:); + dd = tmp2(find(tmp2)); + n2 = n*(n+1)/2; + Atmptriu = Atmp(idxtriu,:); + Atmptril = Atmp(idxtril,:); + if (norm(Atmptriu-Atmptril,'fro') > 1e-13) + fprintf('\n warning: constraint matrices not symmetric.'); + fprintf('\n matrices are symmetrized.\n'); + Atmptriu = 0.5*(Atmptriu+Atmptril); + end + Avec{idxblk,1} = spdiags(dd,0,n2,n2)*Atmptriu; + Ctmp = c(rowidx+blknnz(deblkidx(p))+[1:n*n]); + Ctmp = mexmat(pblk,Ctmp,1); + C{idxblk,1} = 0.5*(Ctmp+Ctmp'); + perm{idxblk,1} = deblkidx(p); + end + end + spblkidx = find(blksize <= smallblkdim); + if ~isempty(spblkidx) + cnt = 0; cnt2 = 0; + spblksize = blksize(spblkidx); + nn = sum(spblksize.*spblksize); + nn2 = sum(spblksize.*(spblksize+1)/2); + pos = zeros(nn,1); + dd = zeros(nn2,1); + idxtriu = zeros(nn2,1); + idxtril = zeros(nn2,1); + for p = 1:length(spblkidx) + n = blksize(spblkidx(p)); + n2 = n*(n+1)/2; + pos(cnt+[1:n*n]) = rowidx+blknnz(spblkidx(p))+[1:n*n]; + %% + %% column-wise positions of upper triangular part + %% + tmp = triu(ones(n)); tmp = tmp(:); + idxtriu(cnt2+[1:n2]) = cnt+find(tmp); + %% + %% row-wise positions of lower triangular part + %% + tmp = tril(reshape([1:n*n],n,n)); tmp = tmp(:); + idxtmp = find(tmp); + [dummy,idxsub] = sort(rem(tmp(idxtmp),n)); + idxtril(cnt2+[1:n2]) = cnt+[idxtmp(idxsub(n+1:end));idxtmp(idxsub(1:n))]; + %% + tmp2 = sqrt(2)*triu(ones(n),1) + speye(n,n); + tmp2 = tmp2(:); + dd(cnt2+[1:n2]) = tmp2(find(tmp2)); + cnt = cnt + n*n; + cnt2 = cnt2 + n2; + end + idxblk = idxblk + 1; + blk{idxblk,1} = 's'; blk{idxblk,2} = blksize(spblkidx); + Atmp = At(pos,:); + Atmptriu = Atmp(idxtriu,:); + Atmptril = Atmp(idxtril,:); + if (norm(Atmptriu-Atmptril,'fro') > 1e-13) + fprintf('\n warning: constraint matrices not symmetric.'); + fprintf('\n matrices are symmetrized.\n'); + Atmptriu = 0.5*(Atmptriu+Atmptril); + end + Avec{idxblk,1} = spdiags(dd,0,length(dd),length(dd))*Atmptriu; + Ctmp = c(pos); + Ctmp = mexmat(blk(idxblk,:),Ctmp,1); + C{idxblk,1} = 0.5*(Ctmp+Ctmp'); + perm{idxblk,1} = spblkidx; + end + end +%% +%%******************************************************************* diff --git a/spotopt/tests/example_Chebyshev.m b/spotopt/tests/example_Chebyshev.m new file mode 100644 index 0000000..df47e36 --- /dev/null +++ b/spotopt/tests/example_Chebyshev.m @@ -0,0 +1,18 @@ +function [pr,objective,answer] = example_Chebyshev() + N=101; + d = 20; + xx = linspace(-1,1,N).'; + Phi = repmat(xx,1,d+1).^repmat((0:d),N,1); + + pr = spotprog; + [pr,c] = pr.newFree(d+1); + + err = abs(xx) - Phi*c; + + [pr,t] = pr.newFree(1); + + pr = pr.withPos(t - err); + pr = pr.withPos(t + err); + answer = 0.0141; + objective = t; +end \ No newline at end of file diff --git a/spotopt/tests/example_RLor2D.m b/spotopt/tests/example_RLor2D.m new file mode 100644 index 0000000..d1eb370 --- /dev/null +++ b/spotopt/tests/example_RLor2D.m @@ -0,0 +1,18 @@ +function [pr,objective,answer] = example_RLor2D() + N = 10; + answer = 0; + pr = spotprog; + [pr,f] = pr.newFree(N); + for i = 1:N + A = randn(2,2); + A = A + A'; + [V,D] = eig(A); + answer = answer + norm(A-V*(D.*(D>=0))*inv(V),'fro'); + + [pr,r] = pr.newRLor(3); + P = [ sqrt(2)*r(1) r(3) ; r(3) sqrt(2)*r(2) ]; + pr = pr.withLor([f(i);A(:)-P(:)]); + end + objective = 1+sum(f); + answer = answer + 1; +end \ No newline at end of file diff --git a/spotopt/tests/example_maxEig.m b/spotopt/tests/example_maxEig.m new file mode 100644 index 0000000..d8f40f6 --- /dev/null +++ b/spotopt/tests/example_maxEig.m @@ -0,0 +1,10 @@ +function [pr,objective,answer] = example_maxEig(n) + A = randn(n,n); + P = A'*A; + answer = max(eig(P)); + + pr = spotprog; + [pr,tau] = pr.newFree(1); + pr = pr.withPSD(tau*eye(n)-P); + objective = tau; +end \ No newline at end of file diff --git a/spotopt/tests/example_sdpProjection.m b/spotopt/tests/example_sdpProjection.m new file mode 100644 index 0000000..70dd33f --- /dev/null +++ b/spotopt/tests/example_sdpProjection.m @@ -0,0 +1,18 @@ +function [pr,objective,answer] = example_sdpProjection(n) + A = randn(n,n); + A = A + A'; + [V,D] = eig(A); + + answer = norm(A-V*(D.*(D>=0))*inv(V),'fro'); + + pr = spotprog; + [pr,p] = pr.newFree(nchoosek(n+1,2)); + P = mss_v2s(p); + pr = pr.withPSD(P); + + [pr,l0] = pr.newFree(1); + + pr = pr.withLor([l0 ; A(:)-P(:)]); + + objective = l0; +end \ No newline at end of file diff --git a/spotopt/tests/example_sosExampleOne.m b/spotopt/tests/example_sosExampleOne.m new file mode 100644 index 0000000..7255997 --- /dev/null +++ b/spotopt/tests/example_sosExampleOne.m @@ -0,0 +1,24 @@ +% Build a homogeneous polynomial, then test the sphere. +n = 4; +x = msspoly('x',n); + +d = 10; +mn = monomials(x,d); + +p = rand(length(mn),1)'*mn; + + +pr = spotsosprog; +[pr,f] = pr.newFree(1); +[pr,q] = pr.newFreePoly(monomials(x,0:d-2)); +pr = pr.withSOS(p + f + q*(1-x'*x)); + +solver = spotsolversedumi; + +sol = pr.minimize(solver,f); + +N = 1000; +X = randn(n,N); +X = X./(repmat(sqrt(sum(X.^2,1)),n,1)); + +[min(msubs(p,x,X)) -sol.eval(f)] \ No newline at end of file diff --git a/spotopt/tests/runExamples.m b/spotopt/tests/runExamples.m new file mode 100644 index 0000000..5308744 --- /dev/null +++ b/spotopt/tests/runExamples.m @@ -0,0 +1,24 @@ +function [] = runExamples(solver) + if nargin < 1, + solver = spotsolversedumi(); + end + + [err,answer,pinf,dinf] = runFeasibleExample(solver,... + @example_maxEig,20); + + if abs(err)/abs(answer) > 1e-6, error('Test failed.'); end + + + [err,answer,pinf,dinf] = runFeasibleExample(solver,... + @example_sdpProjection,20); + if (abs(err)/abs(answer)) > 1e-7, error('Test failed.'); end + + [err,answer,pinf,dinf] = runFeasibleExample(solver,... + @example_Chebyshev); + if (abs(err)/abs(answer)) > 2e-2, error('Test failed.'); end + + [err,answer,pinf,dinf] = runFeasibleExample(solver,... + @example_RLor2D); + if (abs(err)/abs(answer)) > 1e-8, error('Test failed.'); end + +end \ No newline at end of file diff --git a/spotopt/tests/runFeasibleExample.m b/spotopt/tests/runFeasibleExample.m new file mode 100644 index 0000000..22a1a8c --- /dev/null +++ b/spotopt/tests/runFeasibleExample.m @@ -0,0 +1,14 @@ +function [err,answer,pinf,dinf] = runFeasibleExample(solver,test,varargin) + [pr,objective,answer] = test(varargin{:}); + + sol = pr.minimize(solver,objective); + + pinf = sol.primalInfeasible; + dinf = sol.dualInfeasible; + + if pinf || dinf + err = Inf; + else + err = answer-double(sol.eval(objective)); + end +end diff --git a/spotopt/tests/runFeasibleTest.m b/spotopt/tests/runFeasibleTest.m new file mode 100644 index 0000000..99433cf --- /dev/null +++ b/spotopt/tests/runFeasibleTest.m @@ -0,0 +1,18 @@ +function [err,answer,pinf,dinf,derr] = runFeasibleTest(solver,test,varargin) + [pr,objective,answer] = test(varargin{:}); + + [dl,dobj] = toDual(pr,objective); + + sol = solver.minimize(pr,objective); + dsol = solver.minimize(dl,-dobj); + + pinf = sol.primalInfeasible; + dinf = sol.dualInfeasible; + + if pinf || dinf + err = Inf; + else + err = [answer-double(sol.eval(objective)) ... + answer-double(dsol.eval(dobj))]; + end +end diff --git a/spotopt/tests/runTests.m b/spotopt/tests/runTests.m new file mode 100644 index 0000000..50f45d8 --- /dev/null +++ b/spotopt/tests/runTests.m @@ -0,0 +1,15 @@ +function [] = runTests(solver) + if nargin < 1, + solver = spotsolversedumi(); + end + + [err,answer,pinf,dinf] = runFeasibleTest(solver,... + @test_maxEig,spotsdp,20); + + if any(abs(err)/abs(answer) > 1e-6), error('Test failed.'); end + + + [err,answer,pinf,dinf] = runFeasibleTest(solver,... + @test_sdpProjection,spotsdp,20); + if any((abs(err)/abs(answer)) > 1e-7), error('Test failed.'); end +end \ No newline at end of file diff --git a/spotopt/tests/sdp_test.m b/spotopt/tests/sdp_test.m new file mode 100644 index 0000000..b5c414f --- /dev/null +++ b/spotopt/tests/sdp_test.m @@ -0,0 +1,21 @@ +pr = spotsdp; + + +n = 10; +m = 4; +A = randn(n,m); + +[pr,p] = pr.newFree(1); +[pr] = pr.withPos(p); +[pr] = pr.withPSD([ p*eye(n) A ; A' p*eye(m)]); + +solver = spotsolversedumi(struct('fid',1)); + +sol = solver.minimize(pr,p); +sol = solver.minimizeDualForm(pr,p); + +solver = spotsolversdpnal(); +sol = solver.minimizeDualForm(pr,p); + +solver = spotsolversdpt3_4(); +sol = solver.minimizeDualForm(pr,p); \ No newline at end of file diff --git a/spotopt/tests/test_maxEig.m b/spotopt/tests/test_maxEig.m new file mode 100644 index 0000000..9663404 --- /dev/null +++ b/spotopt/tests/test_maxEig.m @@ -0,0 +1,9 @@ +function [pr,objective,answer] = test_maxEig(pr,n) + A = randn(n,n); + P = A'*A; + answer = max(eig(P)); + + [pr,tau] = pr.newFree(1); + pr = pr.withPSD(tau*eye(n)-P); + objective = tau; +end \ No newline at end of file diff --git a/spotopt/tests/test_sdpProjection.m b/spotopt/tests/test_sdpProjection.m new file mode 100644 index 0000000..4687ebf --- /dev/null +++ b/spotopt/tests/test_sdpProjection.m @@ -0,0 +1,17 @@ +function [pr,objective,answer] = test_sdpProjection(pr,n) + A = randn(n,n); + A = A + A'; + [V,D] = eig(A); + + answer = norm(A-V*(D.*(D>=0))*inv(V),'fro'); + + [pr,p] = pr.newFree(nchoosek(n+1,2)); + P = mss_v2s(p); + pr = pr.withPSD(P); + + [pr,l0] = pr.newFree(1); + + pr = pr.withLor([l0 ; A(:)-P(:)]); + + objective = l0; +end \ No newline at end of file diff --git a/tests/spot_build_gram_basis_test.m b/tests/spot_build_gram_basis_test.m new file mode 100644 index 0000000..0f2ba8e --- /dev/null +++ b/tests/spot_build_gram_basis_test.m @@ -0,0 +1,79 @@ +function testPass = spot_build_gram_basis_test + + testPass = []; + + x = msspoly('x',2); + u = msspoly('u',1); + + f = u(1)*x(1)^4+x(1)^2*x(2)^2-u(1)*x(2)^2; + + [F,g,pow] = DecompPoly(f,u(1)); + + m = spot_build_gram_basis(pow,[],F,g); + + if any(m~=[1,1]) + testPass = [testPass,0]; + spotTestFail('Test failed.') + else + testPass = [testPass,1]; + end + + + f = u(1)*x(1)^4+x(1)^2-u(1)*x(2)^2; + [F,g,pow] = DecompPoly(f,u(1)); + m = spot_build_gram_basis(pow,[],F,g); + + if any(m~=[1,0]) + testPass = [testPass,0]; + spotTestFail('Test failed.') + else + testPass = [testPass,1]; + end + + + u = msspoly('u',2); + f = u(1)*x(1)^4+x(1)^2-u(1)*x(2)^2+u(2)*x(1)^6+u(2)*x(2)^3; + [F,g,pow] = DecompPoly(f,u); + m = spot_build_gram_basis(pow,[],F,g); + + if any(m~=[1,0]) + testPass = [testPass,0]; + spotTestFail('Test failed.') + else + testPass = [testPass,1]; + end + + + u = msspoly('u',2); + f = u(2)*x(2)^3; + [F,g,pow] = DecompPoly(f,u); + m = spot_build_gram_basis(pow,[],F,g); + if ~isempty(m) + testPass = [testPass,0]; + spotTestFail('Test failed.') + else + testPass = [testPass,1]; + end + + u = msspoly('u',2); + x = msspoly('x',3); + + f = (u(1)-u(2))*x(1)^4 + u(1)*x(2)^4-u(1)*x(1)^4*x(2)^2+u(2)*x(3)^4+x(3)^2; + [F,g,pow] = DecompPoly(f,u); + m = spot_build_gram_basis(pow,[],F,g); + if any(m~=[0,0,1]) + testPass = [testPass,0]; + spotTestFail('Test failed.') + else + testPass = [testPass,1]; + end + + +end + + +function spotTestFail(string) + error(string) +end + + diff --git a/util/spot_build_gram_basis.m b/util/spot_build_gram_basis.m new file mode 100644 index 0000000..691b5bd --- /dev/null +++ b/util/spot_build_gram_basis.m @@ -0,0 +1,197 @@ +function [mpow] = spot_build_gram_basis(pow,mpow,F,g) + + pow = full(pow); + + if nargin < 2 || isempty(mpow) + mpow = spot_exponent_bound_polytope(pow); + else + mpow = full(mpow); + end + + numHyperPlanes = 2000; + [mpow] = RandomPrune(pow,mpow,numHyperPlanes); + [mpow] = DiagConsistent(pow,mpow); + + if nargin > 2 + mpow = FacialReduction(pow,mpow,F,g); + end + +end + + +%Remove monomials m if 2*m is outside convhull(pow) using +%2*numHyperPlanes random seperating hyperplanes +function mpow = RandomPrune(pow,mpow,numHyperPlanes) + + M = size(pow,2); + + for i=1:numHyperPlanes + + w = randn(M,1); + thres1 = max(pow*w); + thres2 = min(pow*w); + + y = 2*mpow*w; + + %keep mpow if 2*mpow is in convhull(pow) + mpow = mpow(y <= thres1 & y >= thres2,:); + + %Quit if we've pruned mpow enough. + if (size(mpow,1) < 100) + break; + end + + end + +end + +%Get the "diagonally" consistent monomials. 2*mpow must be in pow +%or 2*mpow must be cancelled by a cross term +function [mpow] = DiagConsistent(pow,mpow) + + if (isempty(mpow)) + return + end + + %sorting and casting seems to speed up intersect() + dataType = 'int16'; + pow = sortrows(cast(pow,dataType)); + mpow = sortrows(cast(mpow,dataType)); + + while (1) + + [mext,mextIndx] = ComputeExtremeMonomials(mpow); + [~,mextDeleteIndx] = setdiff(2*mext,pow,'rows'); + + if (~isempty(mextDeleteIndx)) + mpow(mextIndx(mextDeleteIndx),:) = []; + else + break; + end + + end + + mpow = double(mpow); + +end + +%mext: monomials that aren't the midpoints of other monomials +%generalized notion of "extreme point" +function [mext,mextIndx,crossTerms] = ComputeExtremeMonomials(mpow) + + if (size(mpow,1) == 1) + mext = mpow; + mextIndx = 1; + crossTerms = []; + return; + end + + crossTerms = CrossTerms(mpow); + [~,mextIndx] = setdiff(mpow*2,crossTerms,'rows'); + mext = mpow(mextIndx,:); + +end + +%Compute A+A, excluding term A_i + A_j if i = j +function [S] = CrossTerms(A) + + numTerm = size(A,1); + numVar = size(A,2); + + %Get index position of upper triangular part, minus diagonal + posKeep = find(triu(ones(numTerm,numTerm,'uint8'),1)); + + if ~isempty(posKeep) + %Compute A(i,:)+A(j,:) for all i \ne j + S=zeros(length(posKeep),numVar,class(A)); + + for k=1:numVar + temp=bsxfun(@plus,A(:,k),A(:,k)'); + S(:,k) = temp(posKeep); + end + + else + S = []; + end + +end + + +%Apply algorithm from Permenter, Parrilo CDC 2014. +function [mpow] = FacialReduction(pow,mpow,F,g) + + if ~exist('linprog','file') + return + end + + if (isempty(mpow)) + return + end + + powinput = pow; + Finput = F; + ginput = g; + + while(1) + + pow = powinput; + [mext,mextIndx,minkSum] = ComputeExtremeMonomials(mpow); + + %The only relevant monomials are given pow and 2*mext. + %Add 2*mext to pow if necessary and use updated pow to index + %all computations. + powadd = setdiff(2*mext,pow,'rows'); + pow = [powinput;powadd]; + F = [Finput;sparse(size(powadd,1),size(Finput,2))]; + g = [ginput;sparse(size(powadd,1),1)]; + + + [~,facePerpIndx] = setdiff(pow,[mext*2;minkSum],'rows'); + [~,faceDualIndx,faceDualMonomIndx] = intersect(pow,mext*2,'rows'); + + Ffree = F(facePerpIndx,:); + gfree = g(facePerpIndx,:); + Fnneg = F(faceDualIndx,:); + gnneg = g(faceDualIndx,:); + + %Build and solve lp + l = [-Inf*ones(length(facePerpIndx),1);zeros(length(faceDualIndx),1)]; + u = Inf*ones(length(facePerpIndx)+length(faceDualIndx),1); + + sumEq = [zeros(1,length(facePerpIndx)),ones(1,length(faceDualIndx))]; + + A = [[Ffree;Fnneg]';[gfree;gnneg]';sumEq]; + b = [zeros(size(A,1)-1,1);1]; + + [x,~,flag] = linprog(zeros(length(u),1),[],[],A,b,l,u); + + FLAG_SOLUTION_FOUND = 1; + eps = 10^-5; + + %update or quit + if (flag == FLAG_SOLUTION_FOUND) + xdual = x(length(facePerpIndx)+1:end); + deleteIndx = faceDualMonomIndx(xdual > eps); + mpow(mextIndx(deleteIndx),:) = []; + + if (isempty(mpow)) + break + end + + else + break; + end + + end + + +end + + + + + + + + + diff --git a/util/spot_decomp_linear.m b/util/spot_decomp_linear.m new file mode 100644 index 0000000..04d4b6d --- /dev/null +++ b/util/spot_decomp_linear.m @@ -0,0 +1,19 @@ +function [A,b] = spot_decomp_linear(lin,vall) + [veq,peq,Ceq] = decomp(lin); + constant = ~any(peq~=0,2);%all(peq == 0,2); + cnsti = find(constant); + + if isempty(cnsti) + b = sparse(size(Ceq,1),1); + else + b = -Ceq(:,cnsti); + end + + Aeq = Ceq(:,~constant)*peq(~constant,:); + + veqIndices = match(vall,veq); + + % T*vall = veq; + T = sparse(1:length(veq),veqIndices,ones(length(veq),1),length(veq),length(vall)); + A = Aeq*T; +end \ No newline at end of file diff --git a/util/spot_exponent_bound_polytope.m b/util/spot_exponent_bound_polytope.m new file mode 100644 index 0000000..29a364f --- /dev/null +++ b/util/spot_exponent_bound_polytope.m @@ -0,0 +1,18 @@ +function [ptope] = spot_exponent_bound_polytope(pow) +% +% [ptope] = spot_exponent_bound_polytope(pow) +% +% pow -- matrix of nonnegative integers. +% ptope -- matrix of nonnegative integers such that +% ZZ^n intersect conv(pow(:,1),pow(:,2),...,pow(:,end)) +% belongs to ptope. +% + options.sos.newton = 1; + options.sos.csp = 0; + options.verbose = 1; + options.sos.inconsistent = 0; + [C,csclasses] = corrsparsity(pow,options); + mpow = monomialgeneration(pow,csclasses); + + ptope = mpow{1}; +end \ No newline at end of file