Using Prolog & CLP(R) for a system of constraints

Question

I'm looking to use Prolog to generate a random vector that satisfies a system of constraints.

As an example, our user might provide our software with the following information at runtime:

Given a vector <x1, x2, x3, ... x30>, we might have two constraints:

x1 > x2 + x3 + x4
x5 <= sin(x6 + x7)

what I'd like to do is generate a Prolog program that loosely follows the form:

:- random(0.0, 1.0, X1)
:- random(0.0, 1.0, X2)
#...
# its also concievable that the bounds are different than 0 to 1
:- random(0.0, 1.0, X30)

clp(r) :- constraints { 
   X1 > X2 + X3 + X4,
   X5 <= sin(X6 + X7)   
}

?- [ X1, X2, X3, X4, ... X30 ]

which would output a uniformly random vector in the 30-dimensional space.

Is this feasible with Prolog?

There's also the issue of consuming that output. What I'd like to do is have any call to next() re-generate a new vector. Specifically I need to avoid recompilation, since I would like to be able to generate roughly 10,000 of these vectors per second. Can I achieve this level of performance?

I'm hoping to use an embedded (in-process) SWI-Prolog instance on the JVM that the rest of our software runs on. Will that be sufficient?

Answer 1

It's OK!

The approach is OK in principle, and personally, I think Prolog is a good choice for such tasks.

However, there are several subtleties that you need to solve.

First, let us get the syntax of CLP(R) right:

vector([X1,X2,X3,X4,X5,X6,X7]) :-
        { X1 > X2 + X3 + X4,
          X5 =< sin(X6 + X7) }.

Note in particular the use of =< and also the correct use of {}/1 to denote CLP(R) constraints. The token <= is eschewed in Prolog arithmetic because it looks like an arrow that commonly denotes implication in provers.

This is enough to get first answers, even if they are not yet instantiated to concrete solutions:

?- vector(Ls).
Ls = [_1028, _1034, _1040, _1046, _1052, _1058, _1064],
{_1046=_1028-_1034-_1040-_1088, _1088>0.0},
{_1052-sin(_1058+_1064)=<0.0},
{_1052-sin(_1058+_1064)=<0.0},
{_1052-sin(_1058+_1064)=<0.0}.

Using random/1 we can assign random floating point numbers from (0,1) to any of the variables. For example:

?- vector([A,B,C,D,E,F,G]),
   random(A),
   random(B).
A = 0.33797712696696053,
B = 0.7039688010209147,
{D= -0.3659916740539542-C-_894, _894>0.0},
{E-sin(F+G)=<0.0},
{E-sin(F+G)=<0.0},
{E-sin(F+G)=<0.0}.

This solves one part of the task. However, this method fails us in cases like:

?- vector([A,B,C,D,E,F,G]),
   random(A),
   random(B),
   random(C),
   random(D).
false.

Here, the (deterministic!) random number generation conflicts with constraints. There are several ways to get around this. Before I show them, let us constrain the variables to the desired interval, using for example the following definition:

zero_to_one(X) :- { 0 =< X, X =< 1 }.

We can simple state this constraint as one additional requirement:

?- vector([A,B,C,D,E,F,G]),
   maplist(zero_to_one, [A,B,C,D,E,F,G]),
   random(A),
   random(B),
   random(C).

This again yields false.

Method 1: More of the same

One way to solve the above issue is to simply repeat the randomized assignment until a solution is found on backtracking:

?- vector([A,B,C,D,E,F,G]),
   maplist(zero_to_one, [A,B,C,D,E,F,G]),
   random(A),
   random(B),
   repeat,
      random(C).
A = 0.9433451780634803,
B = 0.15859272177823736,
C = 0.706502025956454,
{D>=0.0, _2064=0.07825043032878898-D, D<0.07825043032878898},
{E>=0.0, E=<1.0, F>=0.0, F==0.0, G=<1.0, E-sin(...)=<0.0},
{E-sin(F+G)=<0.0},
{E-sin(F+G)=<0.0},
{E-sin(F+G)=<0.0} .

Thus, we are one step closer to a concrete solution, by which we mean a full instantiation of the vector. The drawbacks are pretty obvious: In extreme cases, we will never find a valid assignment in this way. With slightly better luck, it may take many attempts to find a concrete value for even a single additional variable.

Method 2: Maximize or minimize

A different way to solve this is to use maximize/1 and/or minimize/1 from CLP(R) to use the constraint solver itself to obtain concrete solutions. This only works for linear constraints, and not even for all of those. For example, consider the following queries:

?- { X = sin(Y) },
   maplist(zero_to_one, [X,Y]),
   maximize(X).
false.

And even:

?- { X < 1 }, maximize(X).
false.

Though in contrast:

?- { X =< 1 }, maximize(X).
X = 1.0 .

Now, let us use the following trick to get rid of all nonlinear constraints: We simply assign random floats to X6 and X7, using for example:

?- vector(Ls),
   Ls = [A,B,C,D,E,F,G],
   maplist(zero_to_one, Ls),
   random(F), random(G).

Building on this, we can write:

?- vector(Ls),
   Ls = [A,B,C,D,E,F,G],
   maplist(zero_to_one, Ls),
   random(F), random(G),
   maximize(A), minimize(B+C+D+E).
Ls = [1.0, 0.0, 0.0, 0.0, 0.0, 0.9702069686491169, 0.13220925936558517],
A = 1.0,
B = C, C = D, D = E, E = 0.0,
F = 0.9702069686491169,
G = 0.13220925936558517 .

Thus, we have obtained a concrete solution that satisfies all constraints and has some random components.

Closing remarks

First, to repeat, I think Prolog is a good choice for such tasks. Pruning by the constraint solver can help to eliminate large portions of the search space, and the constraint solver itself can help you to obtain concrete solutions via minimization and maximization. Second, there are still a few issues you need to keep in mind:

• First, solutions generated in this way (by either method) are not random in the sense that each solution is equally likely. Rather, there may be clusters of solutions that are more likely to appear than others.
• As shown above, the equations may require some additional reasoning and experimentation, both to reduce them to linear equations and also to apply the applicable optimization direction. Prolog is well suited for such reasoning, and you can use it to try different strategies easily.
• You may have to find a working trade-off between randomization and deterministic optimization to instantiate remaining vector components. The trade-off may also depend on the entaglement of the vector's components.

Lastly, a very important remark: Implicit random states run counter to the properties we expect from logical relations, since they can cause your predicates to behave quite differently on subsequent invocations, making debugging and systematic testing a nightmare. Therefore, I strongly recommend you make provisions for a random seed, or carry along an explicit state of the random number generator through your code. This will help you to better understand the behaviour of your program and make it completely deterministic. You can later vary the seed to generate different collections of solutions.