Convex Optimization in Java with Project Panama

Project Panama’s Foreign Function and Memory API is the new way of connecting Java programs to non-Java components such as native libraries written in C or C++. This article is about a PoC for calling the convex optimization software ECOS written in C from the Java Virtual Machine with Panama.

Convex optimization is about mathematical optimization of convex functions over convex sets. Many real world problems in a wide range of disciplines such as finance, signal processing, or automatic control systems may be formulated as convex optimization problems.

The PoC has been tested on Windows 10 and Ubuntu 20.04 (WSL2) using Java 17 (Zulu). All code may be found in the PoC’s Github repository.

Call C Code from Java with Panama

In contrast to the Java Native Interface (JNI) which requires to write the native library stubs in C or C++, Panama is a pure Java solution. It allows to access native code via Java method handles. This radically simplifies writing, building, and distributing Java libraries which depend upon native code. At the same time Panama provides performance that is comparable to, sometimes even better than, JNI.

As an example of using Panama, below is a Java unit test for calling the getpid C function from Java

@Test
void panamaSmokeTest() throws Throwable {
    var getpidFuncName = System.getProperty("os.name").startsWith("Windows") ? "_getpid" : "getpid";
    var getpidSymbol = CLinker.systemLookup().lookup(getpidFuncName).orElseThrow();
    var getpidType = MethodType.methodType(int.class);
    var getpidFuncDesc = FunctionDescriptor.of(CLinker.C_INT);
    var getpidHandle = CLinker.getInstance().downcallHandle(getpidSymbol, getpidType, getpidFuncDesc);

    var pid = (int) getpidHandle.invokeExact();

    System.out.println("Process id: " + pid);
    assertTrue(pid > 0);
}

First, CLinker.systemLookup().lookup(getpidFuncName) looks up the symbol with the given name via a system lookup which is suitable to find symbols in the standard C libraries such as getpid and returns the memory address associated with the symbol. Then, CLinker.getInstance().downcallHandle(getpidSymbol, getpidType, getpidFuncDesc) obtains a foreign method handle with the given type and featuring the given function descriptor, which may be used to call the target foreign function at the given address. Finally, (int) getpidHandle.invokeExact() invokes the method handle and casts the return value to an int.

The following command will run the test

"C:\Program Files\Zulu\zulu-17\bin\java.exe" --enable-native-access=ALL-UNNAMED --add-modules jdk.incubator.foreign ...
WARNING: Using incubator modules: jdk.incubator.foreign
Process id: 3864

where --add-modules is required as Panama is still incubating, and --enable-native-access=ALL-UNNAMED allows classpath based code to invoke native code.

Automatic Stubs Generation

In the getpid example, we had to work with C symbol name, system lookup, and method handle in order to call a simple C function. This is a lot of work, but fortunately Panama provides a tool called jextract that greatly simplifies this. jextract automatically generates Java stubs from one or more native C headers.

Currently, the tool may be obtained from one of the Panama Early-Access Builds or directly from the Panama Github repository. We used the Java 17 Panama Early-Access Build.

In order to use jextract for the getpid example, it is enough to create a C header file getpid.h with the following content

int getpid();

Then from the project root simply run jextract with the following command

$JAVA_HOME/bin/jextract \
    -d src/main/java 
    -t com.ustermetrics.getpid.stubs 
    --source 
    getpid.h

where JAVA_HOME should point to the Panama Early-Access Build. This command generates the Java stubs. Additionally, a small manual change in the stubs code is required in order to handle the platform specific symbol name.

The Java unit test for calling getpid simplifies then to

import static com.ustermetrics.getpid.stubs.getpid_h.getpid;

@Test
void jextractSmokeTest() {
    var pid = getpid();

    System.out.println("Process id: " + pid);
    assertTrue(pid > 0);
}

and running the test prints

"C:\Program Files\Zulu\zulu-17\bin\java.exe" --enable-native-access=ALL-UNNAMED --add-modules jdk.incubator.foreign ...
WARNING: Using incubator modules: jdk.incubator.foreign
Process id: 10408

Convex Optimization with ECOS

ECOS is a convex optimization software written in C for solving second-order cone programs (SOCPs). Many convex optimization problems may be formulated as SOCPs, and modern conic optimizers such as ECOS provide a robust and efficient way of solving such problems.

Mathematically, ECOS solves SOCPs in the standard form

\[\begin{align*} & \text{minimize} & & c^T x \\ & \text{subject to} & & A x = b \\ & & & G x + s = h, \qquad s \in \mathcal{K} \end{align*}\]

where \(x \in \mathbf{R}^n\) are the primal variables, \(s \in \mathbf{R}^m\) are slack variables, \(c \in \mathbf{R}^n, \, A \in \mathbf{R}^{p \times n}, \, b \in \mathbf{R}^p, \, G \in \mathbf{R}^{m \times n}, \, h \in \mathbf{R}^m\) are the problem data, and \(\mathcal{K}\) is the convex cone. The cone \(\mathcal{K}\) is the Cartesian product of the following primitive cones

\[\begin{align*} & \text{positive orthant cone} & & \{s \; | \; s \geq 0 \} \\ & \text{second-order cone} & & \{ (t, s) \in \mathbf{R} \times \mathbf{R}^k \; | \; \| s \|_2 \leq t \} \\ & \text{exponential cone} & & \{ (x, y, z) \in \mathbf{R}^3 \; | \; y e^{x/y} \leq z, \, y > 0 \} \end{align*}\]

As an example of convex optimization in quantitative finance consider the Markowitz portfolio optimization problem, where a long-only investor wishes to maximize the expected portfolio return given a limit on the portfolio risk

\[\begin{align*} & \text{maximize} & & \mu^T x \\ & \text{subject to} & & x^T \Sigma x \leq \sigma^2 \\ & & & \mathbf{1} x = 1 \\ & & & x \geq 0 \end{align*}\]

where \(x\) is the unknown vector of portfolio allocations, \(\mu\) is the estimated expected return vector, \(\Sigma\) is the estimated covariance matrix, and \(\sigma\) is the given limit on the portfolio risk.

In order to solve the Markowitz portfolio optimization problem with ECOS, it is necessary to convert it to a SOCP \[\begin{align*} & \text{minimize} & & -\mu^T x \\ & \text{subject to} & & \mathbf{1} x = 1 \\ & & & x \geq 0 \\ & & & \| G^T x \|_2 \leq \sigma \end{align*}\]

where \(G\) is such that \(\Sigma = G G^T\). \(G\) may be found e.g. by the Cholesky decomposition. The two inequalities may be modeled with a positive orthant and a second-order cone, respectively.

Calling ECOS from Java with Panama

In order to call ECOS from Java the three C functions

pwork* ECOS_setup(idxint n, idxint m, idxint p, idxint l, idxint ncones, idxint* q, idxint nex,
                   pfloat* Gpr, idxint* Gjc, idxint* Gir,
                   pfloat* Apr, idxint* Ajc, idxint* Air,
                   pfloat* c, pfloat* h, pfloat* b);
idxint ECOS_solve(pwork* w);
void ECOS_cleanup(pwork* w, idxint keepvars);

need Java stubs. ECOS_setup defines the dimensions and the problem data of the SOCP and allocates work space for the solver. ECOS_solve solves the SOCP and returns the optimality status from the solver. Finally, ECOS_cleanup frees the allocated memory.

jextract automatically generates these stubs with the following command

$JAVA_HOME/bin/jextract \
    -C "-D DLONG" \
    -C "-D LDL_LONG" \
    -C "-DSuiteSparse_long=long long" \
    -d src/main/java \
    -l ecos \
    -t com.ustermetrics.ecos.stubs \
    --source \
    -I $ECOS_HOME/external/SuiteSparse_config \
    --include-function ECOS_setup \
    --include-function ECOS_solve \
    --include-function ECOS_cleanup \
    --include-function ECOS_ver \
    --include-struct pwork \
    --include-struct stats \
    $ECOS_HOME/include/ecos.h

where we added some extra arguments such as pass through arguments for clang and names of struct definitions to include in the stubs.

Using the stubs and Panama, Java code for solving SOCPs with ECOS may be compactly written as

public static Solution solve(double[] c, double[] Apr, long[] Ajc, long[] Air, double[] b, double[] Gpr,
                             long[] Gjc, long[] Gir, double[] h, long l, long[] q, long nex) {
    var n = c.length;
    var m = h.length;
    var p = b.length;
    var ncones = q.length;

    long exitCode;
    double[] solution = null;
    double cost = Double.NaN;

    try (var sc = ResourceScope.newConfinedScope()) {
        var alloc = SegmentAllocator.arenaAllocator(sc);

        var qSeg = alloc.allocateArray(C_LONG_LONG, q);
        var GprSeg = alloc.allocateArray(C_DOUBLE, Gpr);
        var GjcSeg = alloc.allocateArray(C_LONG_LONG, Gjc);
        var GirSeg = alloc.allocateArray(C_LONG_LONG, Gir);
        var AprSeg = alloc.allocateArray(C_DOUBLE, Apr);
        var AjcSeg = alloc.allocateArray(C_LONG_LONG, Ajc);
        var AirSeg = alloc.allocateArray(C_LONG_LONG, Air);
        var cSeg = alloc.allocateArray(C_DOUBLE, c);
        var hSeg = alloc.allocateArray(C_DOUBLE, h);
        var bSeg = alloc.allocateArray(C_DOUBLE, b);

        var workAddr = ECOS_setup(n, m, p, l, ncones, qSeg, nex, GprSeg, GjcSeg, GirSeg, AprSeg, AjcSeg, AirSeg,
                cSeg, hSeg, bSeg);
        if (NULL.equals(workAddr)) {
            throw new IllegalArgumentException("Something went wrong in ECOS_setup()");
        }

        exitCode = ECOS_solve(workAddr);
        if (exitCode == 0) {
            var workSeg = workAddr.asSegment(pwork.sizeof(), sc);
            var xAddr = pwork.x$get(workSeg);
            var xSeg = xAddr.asSegment(C_DOUBLE.byteSize() * n, sc);
            solution = xSeg.toDoubleArray();

            var infoAddr = pwork.info$get(workSeg);
            var infoSeg = infoAddr.asSegment(stats.sizeof(), sc);
            cost = stats.pcost$get(infoSeg);

            ECOS_cleanup(workAddr, 0);
        }
    }

    return new Solution(exitCode, solution, cost);
}

First, try (var sc = ResourceScope.newConfinedScope()) { ... } creates a resource scope that is used with the try-with-resources construct. The resource scope controls the life-cycle of the memory segments inside the try-with-resources. When a resource scope is closed, access to resources managed by that scope is no longer allowed. Then, with multiple alloc.allocateArray( ... ) blocks of memory that hold the SOCP problem data are allocated and initialized. After calling ECOS_setup and ECOS_solve, the solution to the problem is extracted with

var workSeg = workAddr.asSegment(pwork.sizeof(), sc);
var xAddr = pwork.x$get(workSeg);
var xSeg = xAddr.asSegment(C_DOUBLE.byteSize() * n, sc);
solution = xSeg.toDoubleArray();

Finally, ECOS_cleanup frees the allocated native memory on the C side.

As an example of using the above solver method, below is a Java unit test that uses the dense and sparse vector and matrix classes of Apache Spark’s machine learning library to build, solve and verify the Markowitz portfolio optimization problem as explained above

@Test
void ecosSolveShouldReturnOptimalSolution() {
    var c = new DenseVector(new double[]{-0.05, -0.06, -0.08, -0.06, 0.});
    var A = new DenseMatrix(1, 5, new double[]{1., 1., 1., 1., 0.});
    var b = new DenseVector(new double[]{1.});
    var G = new DenseMatrix(10, 5, new double[]{
            -1., 0., 0., 0., 0.,
            0., -1., 0., 0., 0.,
            0., 0., -1., 0., 0.,
            0., 0., 0., -1., 0.,
            0., 0., 0., 0., 1.,
            0., 0., 0., 0., -1.,
            -0.15, -0.02, -0.1, -0.15, 0.,
            0., -0.198997487421324, -0.16583123951776996, -0.10552897060221729, 0.,
            0., 0., -0.158113883008419, -0.17392527130926083, 0.,
            0., 0., 0., -0.16159714218895202, 0.
    }, true);
    var h = new DenseVector(new double[]{0., 0., 0., 0., 0.2, 0., 0., 0., 0., 0.});
    var l = 5;
    var q = new long[]{5};
    var nex = 0;

    var As = A.toSparseColMajor();
    var Gs = G.toSparseColMajor();

    var solution = EcosSolver.solve(c.values(), As.values(), toLongArray(As.colPtrs()),
            toLongArray(As.rowIndices()), b.values(), Gs.values(), toLongArray(Gs.colPtrs()),
            toLongArray(Gs.rowIndices()), h.values(), l, q, nex);

    System.out.println(solution);
    assertEquals(0, solution.exitCode());
    var tol = Math.ulp(1.0); // Machine epsilon
    assertArrayEquals(new double[]{0.24879020572078372, 0.049684806182020855, 0.7015249845663684,
            3.5308169265756875e-09, 0.19999999978141014}, solution.solution(), tol);
    assertEquals(-0.07154259763411892, solution.cost(), tol);
}

with

\[ \mu=(0.05, 0.06, 0.08, 0.06)^T \]

\[ \Sigma = \left[\begin{array} {rrr} 0.02250 & 0.00300 & 0.01500 & 0.02250 \\ 0.00300 & 0.04000 & 0.03500 & 0.02400 \\ 0.01500 & 0.03500 & 0.06250 & 0.06000 \\ 0.02250 & 0.02400 & 0.06000 & 0.09000 \end{array}\right] \]

and \(\sigma = 0.2\).

For the second-order cone the Cholesky decomposition of \(\Sigma\) is used, and the cone is modeled as \(\| G^T x \|_2 \leq t, \; t \leq \sigma\), where \(t\) is an additional primal variable. Hence, ECOS solves a five-dimensional problem, where the first four elements of the solution represent the optimal portfolio allocations, and the last element of the solution is \(t\).

Running the test prints

"C:\Program Files\Zulu\zulu-17\bin\java.exe" --enable-native-access=ALL-UNNAMED --add-modules jdk.incubator.foreign ...
WARNING: Using incubator modules: jdk.incubator.foreign
Solution{exitCode=0, solution=[0.24879020572078372, 0.049684806182020855, 0.7015249845663684, 3.5308169265756875E-9, 0.19999999978141014], cost=-0.07154259763411892}

ECOS 2.0.10 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  -6.235e-02  -2.658e-01  +7e+00  6e-01  6e-01  1e+00  1e+00    ---    ---    1  1  - |  -  - 
 1  -6.268e-02  -1.670e-01  +2e+00  6e-02  7e-02  4e-03  3e-01  0.8551  8e-02   1  1  1 |  0  0
 2  -6.336e-02  -6.999e-02  +2e-01  4e-03  6e-03  1e-03  3e-02  0.9203  2e-02   1  1  1 |  0  0
 3  -7.130e-02  -7.313e-02  +5e-02  9e-04  2e-03  6e-04  8e-03  0.8375  2e-01   1  1  1 |  0  0
 4  -7.140e-02  -7.165e-02  +7e-03  1e-04  2e-04  8e-05  1e-03  0.8980  4e-02   1  1  1 |  0  0
 5  -7.153e-02  -7.158e-02  +1e-03  2e-05  5e-05  2e-05  2e-04  0.8330  5e-02   1  1  1 |  0  0
 6  -7.154e-02  -7.154e-02  +9e-05  2e-06  3e-06  1e-06  1e-05  0.9890  5e-02   1  1  1 |  0  0
 7  -7.154e-02  -7.154e-02  +4e-06  7e-08  1e-07  6e-08  6e-07  0.9592  2e-04   1  1  1 |  0  0
 8  -7.154e-02  -7.154e-02  +1e-07  2e-09  4e-09  2e-09  2e-08  0.9786  1e-02   1  1  1 |  0  0
 9  -7.154e-02  -7.154e-02  +5e-09  9e-11  2e-10  8e-11  8e-10  0.9620  2e-04   1  1  1 |  0  0

OPTIMAL (within feastol=1.7e-10, reltol=6.7e-08, abstol=4.8e-09).
Runtime: 0.000323 seconds.

Summary and Conclusion

This post shows how to use Panama’s Foreign Function and Memory API and its automatic stubs generation tool jextract to call the convex optimization software ECOS written in C from the Java Virtual Machine. As an example, it is shown how to solve the famous Markowitz portfolio optimization problem from quantitative finance in Java.

In contrast to the old Java Native Interface JNI, Panama, in particular in combination with the jextract tool, radically simplifies writing, building, and distributing Java libraries which depend upon native code.

For Java 19 that is scheduled for September 2022 and now under active development, Panama’s Foreign Function and Memory API will transition from an incubating API in module jdk.incubator.foreign to a preview API in java.base. Finally, a list of other examples using Panama and jextract may be found here.

Adrian Trapletti
Adrian Trapletti
PhD, CEO

Quant, software engineer, and consultant mostly investment industry. Long-term contributor and package author R Project for Statistical Computing.