Converting a quadratic program to a conic program (second order cone program)

September 28th, 2014

Mosek, a commercial solver, writes in its documentation that it’s often advantageous to convert a separable quadratic program:

min_x 1/2 |Fx - c|^2
subject to Ax <= b

with x a n-long vector, into a conic program of the form

min_{x,t,y} y - x'F'c
subject to  Ax <= b
and         Fx - t = 0
and         2y >= t't

introducing t as a new vector variable and y as a new scalar variable.

It’s rather straightforward to verify that these two problems are equivalent. But what’s the intuition behind it?

To understand this, we’ll need to be comfortable with two ways of visualizing an energy function: as a line plot or height-field surface and viewed from above as a topographic map of isolines.

Recall that an energy function E(x) is just a scalar function of one or more variables in the vector x. To be sure, this means E(x) takes in a list of numbers in x and spits out a single number.

In our original quadratic program, the energy was a quadratic function of x. If x were a only a single variable (scalar) then we can rewrite 1/2 |Fx - c|^2 as f^2/2 x^2 - 2fc x + c^2. This is simply a quadratic function x, familiar after substitution as: ax^2 + bx + c. We know that plotting this will reveal a parabola:

quadratic energy parabola

Minimizing E(x) in this picture means walking along values of x and tracing the curve until reaching the vertex, the point with minimal E(x) and, notably, zero derivative.

Conceptually when we convert this simpler 1D quadratic problem E(x) = ax^2 + bx + c into a conic program we introduce a new variable y. It’s important that we really recognize y as a variable. So far we don’t know its value. If we stopped here we’d have an optimization problem:

min_{x,y} ax^2 + bx + c

where y spans an infinite null space of solutions: the energy only measures the effect of x; any choice of y is OK.

But here comes the trick, we place y directly in the objective function: E_C(x,y) = y. This new energy is linear in only y. If we know plot it as a topographical map over the xy-plane we see a function decreasing indefinitely in the downward direction:

linear energy topographical map

Obviously if we stop here and minimize over x and y, we could get any value of x and we’ll simply get y=-∞. The final step is to add a conic constraint on y according to the original quadratic function of x. Namely we ask that y >= ax^2 + bx + c: y should be above the original quadratic function. Since we’re only working with scalar x and scalar y we can plot the feasible region as on the xy-plane. y must be above the red region:

conic constraint feasible region

Now we can imagine minimizing y by traveling downward, we’ll eventually hit the red region’s wall. Since we have full freedom to move in the x direction, too, we can continue downward only if we slide along the wall toward the minimum of the original quadratic function E(x).

Higher dimensions: This scenario is significantly easier to visualize for the 1d problem. However, the intuition extends to higher dimensions. If x is a vector then E(x) is a hyper-paraboloid over x. When we introduce y, we still just add a single new orthogonal dimension. Then minimizing the value of y we’ll smash into this hyper-paraboloid and slide toward the minimum.

Solvers: This explanation is meant to give some intuition as to why this conversion works. It’s not meant as a explanation of how interior-point conic programming solvers work. I also do not explain why it’s sometimes more efficient to convert least squares quadratic programs to conic programs. Admittedly the conditions for a successful conversion in terms of efficiency still allude me.

“Skinning Cubic Bezier Splines and Catmull Clark Subdivision Surfaces” preprint available

September 28th, 2014

skinning vector graphics in 2D and 3D

I’ve uploaded a preprint of a paper Songrun Liu, Yotam Gingold and I have recently finished. The paper’s entitled “Skinning Cubic Bezier Splines and Catmull Clark Subdivision Surfaces” and will be published in the upcoming proceedings of SIGGRAPH Asia 2014.

Abstract: Smooth space deformation has become a vital tool for the animation and design of 2D and 3D shapes. Linear methods, under the umbrella term of “linear blend skinning”, are the de facto standard for 3D animations. Unfortunately such approaches do not trivially extend to deforming vector graphics, such as the cubic Bézier splines prevalent in 2D or subdivision surfaces in 3D. We propose a variational approach to reposition the control points of cubic Bézier splines and Catmull-Clark subdivision surfaces—or any linear subdivision curves or surfaces—to produce curves or surfaces which match a linear blend skinning deformation as closely as possible. Exploiting the linearity of linear blend skinning, we show how this optimization collapses neatly into the repeated multiplication of a matrix per handle. We support C0, C1, G1, and fixed-angle continuity constraints between adjacent Bézier curves in a spline. Complexity scales linearly with respect to the number of input curves and run-time performance is fast enough for real-time editing and animation of high-resolution shapes.

Update: We now have a project page, too.

Accompanying video for SIGGRAPH Asia paper “Skinning Cubic Bezier Splines and Catmull Clark Subdivision Surfaces “

September 28th, 2014

Here’s the accompanying video for our upcoming SIGGRAPH Asia paper “Skinning Cubic Bezier Splines and Catmull Clark Subdivision Surfaces” by Songrun Liu, Alec Jacobson and Yotam Gingold.

Also, check out the paper.

2014 SciVis paper “Fast and Memory-Efficient Topological Denoising of 2D and 3D Scalar Fields”

September 22nd, 2014

David Günther, Jan Reininghaus, Hans-Peter Seidel, Olga Sorkine-Hornung, Tino Weinkauf and I will publish our paper
“Fast and Memory-Efficient Topological Denoising of 2D and 3D Scalar Fields” in TVCG and David will present the work at SciVis this year. This paper is a much-needed extension of our earlier SGP paper on optimizing over monotonic functions. This new paper provides two new heuristics: 1) a plaided domain decomposition enables efficient optimization on much larger datasets (notably 3D images, e.g. voxel grids), and 2) iterating on our convexification of the monotonicity constraint. This iteration idea greatly improves the quality of solutions when the gradient direction of the harmonic initial guess is far from matching the optimal solution, particularly important in the case of data smoothing. In some sense, one could view this as a sort of trust-region-ish algorithm, but our problem specific knowledge allows us to make much better progress.

Download the paper.

Using embree compiled with gcc from a program compiled with clang

September 12th, 2014

I compiled Embree using my standard compiler gcc 4.7 (I’m holding out on clang for openmp support), but for a specific project I need to use clang. When I try to link against libembree.a and libsys.a I get this sort of error:

  "std::string::find(char, unsigned long) const", referenced from:
      embree::AccelRegistry::create(std::string, embree::RTCGeometry*) in libembree.a(registry_accel.o)
      embree::BuilderRegistry::build(std::string, embree::TaskScheduler::Event*, embree::Accel*) in libembree.a(registry_builder.o)
      embree::IntersectorRegistry<embree::Intersector1>::get(std::string, std::string, embree::Accel*) in libembree.a(registry_intersector.o)
      embree::IntersectorRegistry<embree::Intersector4>::get(std::string, std::string, embree::Accel*) in libembree.a(registry_intersector.o)
  "std::string::compare(char const*) const", referenced from:
      embree::AccelRegistry::create(std::string, embree::RTCGeometry*) in libembree.a(registry_accel.o)
      embree::BuilderRegistry::build(std::string, embree::TaskScheduler::Event*, embree::Accel*) in libembree.a(registry_builder.o)
      embree::IntersectorRegistry<embree::Intersector1>::get(std::string, std::string, embree::Accel*) in libembree.a(registry_intersector.o)
      embree::IntersectorRegistry<embree::Intersector4>::get(std::string, std::string, embree::Accel*) in libembree.a(registry_intersector.o)
...

I tried using the flag -stdlib=libc++ but this doesn’t fix the problem. Using -stdlib=libstdc++ most errors go away except for:

  "std::ctype<char>::_M_widen_init() const", referenced from:
      embree::rtcBuildAccel(embree::RTCGeometry*, char const*) in libembree.a(embree.o)
      embree::BVHStatisticsT<embree::BVH2>::print() in libembree.a(bvh2.o)
      embree::BVHStatisticsT<embree::BVH4>::print() in libembree.a(bvh4.o)
      embree::BVHBuilderT<embree::BVH4, embree::HeuristicBinning<0> >::createLeafNode(unsigned long, embree::atomic_set<embree::PrimRefBlock>&, embree::HeuristicBinning<0>::PrimInfo const&) in libembree.a(bvh4.o)
      embree::BVHBuilderT<embree::BVH4, embree::HeuristicBinning<2> >::createLeafNode(unsigned long, embree::atomic_set<embree::PrimRefBlock>&, embree::HeuristicBinning<2>::PrimInfo const&) in libembree.a(bvh4.o)
      embree::BVHBuilderT<embree::BVH4, embree::HeuristicSpatial<0> >::createLeafNode(unsigned long, embree::atomic_set<embree::PrimRefBlock>&, embree::HeuristicSpatial<0>::PrimInfo const&) in libembree.a(bvh4.o)
      embree::BVHBuilderT<embree::BVH4, embree::HeuristicSpatial<2> >::createLeafNode(unsigned long, embree::atomic_set<embree::PrimRefBlock>&, embree::HeuristicSpatial<2>::PrimInfo const&) in libembree.a(bvh4.o)

This error seems trickier and to remove it I explicitly link against libgcc_s.a found in my gcc’s libraries: add the linker argument /opt/local/lib/gcc47/libstdc++.a. This gets me close, but there’s a new error:

  "___emutls_get_address", referenced from:
      ___cxa_get_globals_fast in libstdc++.a(eh_globals.o)
      ___cxa_get_globals in libstdc++.a(eh_globals.o)

Turns out that I need another library from gcc making the arguments now: /opt/local/lib/gcc47/libstdc++.a /opt/local/lib/gcc47/libgcc_s.1.dylib. And this finally works.

Update: Oddly it seems on the static libstdc++.a has the final problem. Linking to the dymnamic library libstdc++.dylib in the same directory works fine (maybe because I’m just avoiding the issue with this example).

Eigen + C++11’s auto gotcha

September 10th, 2014

Today marks the second time running into this gotcha involving Eigen and C++11’s auto. It may even be a bug.

Here’s the little example program:

Eigen::MatrixXd X(2,2);                                                    
X<<1,2,
   3,4;                                                                    
const auto mid = (X.colwise().maxCoeff()+X.colwise().minCoeff()).eval()*0.5;  
cout<<mid<<endl;                                                           
X.setZero();                                                               
cout<<mid<<endl;
cout<<(X.colwise().maxCoeff()+X.colwise().minCoeff()).eval()*0.5<<endl;    

First I initialize a little matrix and then take the average of it’s column-wise max and min vectors and store that in mid. I’m being careful to use .eval() because this should truly evaluate the expression. Then I print mid, zap X to zero, print mid again and—to double check—print the expression applied to the new all-zero X:

With clang this produces:

2 3
  1 1.5
0 0

With gcc I get:

2 3
-6.44115e-232  1.34339e+154
0 0

So it seems that after X is set to zeros, mid contains some random bad value.

Changing auto above to char and reading the compiler error message I can reveal the true type of the expression

const ScalarMultipleReturnType {aka const Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, const Eigen::Matrix<double, 1, -1> >}

Okay, so definitely not something simple like Eigen::RowVector<...> which I’d expected. Not sure why this is returning junk when the values of X change. For now, and for this simple example, the easy solution is to explicitly cast to the type I want:

const RowVector2d mid = (X.colwise().maxCoeff()+X.colwise().minCoeff()).eval()*0.5;

Compiling C++ with C++11 support on bluehost servers

September 8th, 2014

I asked bluehost to install a modern version of gcc with C++11 support. They told me that only gcc4.3 was deemed stable.

So then I tried to compile gcc4.7 and gcc4.9 from scratch. This seemed promising but eventually I hit some apparent bug in the make routine of gcc giving me the error:

"cp: cannot stat `libgcc_s.so.1': No such file or directory"

Manually copying the file for make did not help.

So then I switched to trying to compile clang. This worked. I followed the instructions in this answer.

Namely, I added this to my .bashrc

export LD_LIBRARY_PATH=/usr/lib/gcc/x86_64-redhat-linux/4.4.7/:$LD_LIBRARY_PATH
export CC=/usr/bin/gcc
export CXX=/usr/bin/g++

Then issued the following:

source .bashrc
wget http://llvm.org/releases/3.3/cfe-3.3.src.tar.gz
tar xzf llvm-3.3.src.tar.gz && cd llvm-3.3.src/tools/ && tar xzf ../../cfe-3.3.src.tar.gz
cd llvm-3.3.src
mv tools/cfe-3.3.src tools/clang
./configure --prefix=$HOME/llvm-3.3.src/llvm
make -j8
make install

Now to test it out, I created a file test.cpp with this inside:

#include <iostream>
int main(int argc, char *argv[])
{
  const auto & hello = []()
  {
    std::cout<<"Hello, world."<<std::endl;
  };
  hello();
}

Then I could try to compile with:

llvm-3.3.src/llvm/bin/clang++ -std=c++11 test.cpp -o test

But I get the error:

In file included from test.cpp:1:
In file included from /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/iostream:40:
In file included from /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/ostream:40:
In file included from /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/ios:40:
In file included from /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/exception:148:
/usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/exception_ptr.h:143:13: error: unknown type name
      'type_info'
      const type_info*
            ^
1 error generated.

This seems to be a known bug with old “stable” versions of gcc’s standard library. The fix is to add a class declaration before the #include <iostream>:

#ifdef __clang__
    class type_info;
#endif
#include <iostream>
int main(int argc, char *argv[])
{
  const auto & hello = []()
  {
    std::cout<<"Hello, world."<<std::endl;
  };
  hello();
}

This compiles and runs.

Compile Embree without GLUT

September 8th, 2014

Trying to compile embree on a server, I was notified that I didn’t have GLUT installed:

CMake Error at /home2/alecjaco/share/cmake-2.8/Modules/FindPackageHandleStandardArgs.cmake:97 (MESSAGE):
  Could NOT find GLUT (missing: GLUT_glut_LIBRARY GLUT_INCLUDE_DIR)

I shouldn’t need this because I don’t care about visualizing the results. To get around this requirement I simply comment out the following line in the top level CMakeLists.txt:

ADD_SUBDIRECTORY(examples)

becomes

#ADD_SUBDIRECTORY(examples)

CGAL function pointer callback in isosurface example

September 4th, 2014

I really struggled with the CGAL templating madness in their isosurface example today. The following example contours the 0-level isosurface given an implicit function of a sphere (i.e. a pointer to a function sphere_function):

#include <CGAL/Surface_mesh_default_triangulation_3.h>
#include <CGAL/Complex_2_in_triangulation_3.h>
#include <CGAL/make_surface_mesh.h>
#include <CGAL/Implicit_surface_3.h>
// default triangulation for Surface_mesher
typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
// c2t3
typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;
typedef Tr::Geom_traits GT;
typedef GT::Sphere_3 Sphere_3;
typedef GT::Point_3 Point_3;
typedef GT::FT FT;
typedef FT (*Function)(Point_3);
typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;
FT sphere_function (Point_3 p) {
  const FT x2=p.x()*p.x(), y2=p.y()*p.y(), z2=p.z()*p.z();
  return x2+y2+z2-1;
}
int main() {
  Tr tr;            // 3D-Delaunay triangulation
  C2t3 c2t3 (tr);   // 2D-complex in 3D-Delaunay triangulation
  // defining the surface
  Surface_3 surface(sphere_function,             // pointer to function
                    Sphere_3(CGAL::ORIGIN, 2.)); // bounding sphere
  // Note that "2." above is the *squared* radius of the bounding sphere!
  // defining meshing criteria
  CGAL::Surface_mesh_default_criteria_3<Tr> criteria(30.,  // angular bound
                                                     0.1,  // radius bound
                                                     0.1); // distance bound
  // meshing surface
  CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Non_manifold_tag());
  std::cout << "Final number of points: " << tr.number_of_vertices() << "\n";
}

I wanted to adapt this example to use a different implicit function, so I tried creating a sphere function that also took another parameter (x-axis scaling):

FT scale_sphere_function (double a,Point_3 p) {
  const FT x2=a*p.x()*p.x(), y2=p.y()*p.y(), z2=p.z()*p.z();
  return x2+y2+z2-1;
}

Then my idea was to use std::bind to bind a local variable or value to the first argument so I could pass it to Surface_3. Something like

double a = 2.0;
Surface_3 surface(
  std::bind(&scale_sphere_function,a,std::placeholders::_1),
  Sphere_3(CGAL::ORIGIN, 2.)); // bounding sphere

But I get compiler errors because it’s trying to convert the std::bind output (std::function<>) to a function pointer:

/opt/local/include/CGAL/Implicit_surface_3.h:50:5: note: CGAL::Implicit_surface_3<GT, Function_>::Implicit_surface_3(CGAL::Implicit_surface_3<GT, Function_>::Function, CGAL::Implicit_surface_3<GT, Function_>::Sphere_3, CGAL::Implicit_surface_3<GT, Function_>::FT, CGAL::Implicit_surface_3<GT, Function_>::Geom_traits) [with GT = CGAL::Robust_circumcenter_traits_3<CGAL::Epick>; Function_ = double (*)(CGAL::Point_3<CGAL::Epick>); CGAL::Implicit_surface_3<GT, Function_>::Function = double (*)(CGAL::Point_3<CGAL::Epick>); CGAL::Implicit_surface_3<GT, Function_>::Sphere_3 = CGAL::Sphere_3<CGAL::Epick>; CGAL::Implicit_surface_3<GT, Function_>::FT = double; CGAL::Implicit_surface_3<GT, Function_>::Geom_traits = CGAL::Robust_circumcenter_traits_3<CGAL::Epick>]
/opt/local/include/CGAL/Implicit_surface_3.h:50:5: note:   no known conversion for argument 1 from 'std::_Bind_helper<false, double (&)(double, CGAL::Point_3<CGAL::Epick>), int, const std::_Placeholder<1>&>::type {aka std::_Bind<double (*(int, std::_Placeholder<1>))(double, CGAL::Point_3<CGAL::Epick>)>}' to 'CGAL::Implicit_surface_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, double (*)(CGAL::Point_3<CGAL::Epick>)>::Function {aka double (*)(CGAL::Point_3<CGAL::Epick>)}'
/opt/local/include/CGAL/Implicit_surface_3.h:34:9: note: constexpr CGAL::Implicit_surface_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, double (*)(CGAL::Point_3<CGAL::Epick>)>::Implicit_surface_3(const CGAL::Implicit_surface_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, double (*)(CGAL::Point_3<CGAL::Epick>)>&)
/opt/local/include/CGAL/Implicit_surface_3.h:34:9: note:   candidate expects 1 argument, 2 provided
/opt/local/include/CGAL/Implicit_surface_3.h:34:9: note: constexpr CGAL::Implicit_surface_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, double (*)(CGAL::Point_3<CGAL::Epick>)>::Implicit_surface_3(CGAL::Implicit_surface_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, double (*)(CGAL::Point_3<CGAL::Epick>)>&&)
/opt/local/include/CGAL/Implicit_surface_3.h:34:9: note:   candidate expects 1 argument, 2 provided

I tried all sorts of conversions and casts and uses of target() (which is pointless and only works when you originally started with a function pointer) to twist my std::bind into a function pointer. But the problem turned out to be very simple. The example had defined Surface_3 to use a function pointer:

typedef FT (*Function)(Point_3);
typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;

I can just define it to use a std::function<> instead:

typedef std::function<FT (Point_3)> Function;
typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;

Then I can create a Surface_3 with std::bind:

double a = 2.0;
Surface_3 surface(
  std::bind(&sphere_function,a,std::placeholders::_1),
  Sphere_3(CGAL::ORIGIN, 2.)); // bounding sphere

or even an inline lambda with capturing:

Surface_3 surface(
  [&a](Point_3 p)->FT{ return sphere_function(a,p);},
  Sphere_3(CGAL::ORIGIN, 2.)); // bounding sphere

I’m apparently not the only one with problem. Probably CGAL should just use std::function or boost::function in the example since it’s easy to convert to those from a function pointer (but as I learned it’s impossible the other way).

Skinning.org Siggraph 2014 Course Notes

August 12th, 2014

skinning course at siggraph

Zhigang Deng, Ladislav Kavan, JP Lewis and I are teaching a course this week at siggraph.

We’ll cover a wide variety of skinning related techniques and we’re hosting our notes at skinning.org. Check it out and tell us what you think!