Motivating the Template Function System

The template function system is useful for operators that are defined for domains that contain a large number of types.

Consider the function nsolve which implements the proximal operator for the l1 norm.

;; Compute the solution to the following objective
;;
;;   arg min{x}  alpha || x ||_1 + (lambda/2) || m - x ||_2^2
;;
;; where
;;   alpha      a real greater than zero.
;;   lambda     a real greater than zero.
;;   m          a vector of reals.

(defun nsolve (alpha lambda m)
  (declare (type (real (0)) alpha)
           (type (real (0)) lambda)
           (type array m))
  (let* ((element-type (array-element-type m))
         (zero (coerce 0 element-type))
         (pos-k (coerce (/ alpha lambda) element-type))
         (neg-k (- pos-k)))
    (dotimes (i (array-total-size m))
      (let* ((m-i (row-major-aref m i))
             (v-i (cond ((> m-i pos-k)
                         (- m-i pos-k))
                        ((< m-i neg-k)
                         (+ m-i neg-k))
                        (t
                         zero))))
        (setf (row-major-aref m i) v-i))))
  (values))

The arguments alpha and lambda are non negative reals and the argument m is an array of reals.

The function stores the solution x to the proximal operator in the input array m.

The common lisp type real is a super type of a great number of types including the floating point types short-float, single-float, double-float and long-float.

The common lisp specification allows an implementation to provide a unique representation for each floating point type. This will impact the performance of the function nsolve as the addition and subtraction functions must consider the representation at runtime. This is confirmed empirically by evaluating the following code in SBCL.

(defun time-nsolve ()
  (let* ((element-type 'double-float)
         (magnitude (coerce 20 element-type))
         (dimensions 10000)
         (count 3000)
         (alpha 1d0)
         (lambda 1d0)
         (m (make-array dimensions :element-type element-type)))
    (dotimes (i (array-total-size m))
      (let* ((v-i (random magnitude)))
        (setf (row-major-aref m i) (ecase (random 2)
                                     (0 v-i)
                                     (1 (- v-i))))))

    (format t "Timing: ~A" 'nsolve)
    (time
      (dotimes (i count)
        (nsolve alpha lambda m)))))

(time-nsolve)
;;
;; Timing: nsolve
;; Evaluation took:
;;   2.020 seconds of real time
;;   2.016454 seconds of total run time (2.007544 user, 0.008910 system)
;;   [ Run times consist of 0.009 seconds GC time, and 2.008 seconds non-GC time. ]
;;   99.80% CPU
;;   4,040,325,780 processor cycles
;;   709,685,808 bytes consed

To improve the performance of the function nsolve we turn to inlining.

(declaim (inline %nsolve))
(defun %nsolve (alpha lambda m)
  (declare (type (real (0)) alpha)
           (type (real (0)) lambda)
           (type array m))
  (let* ((element-type (array-element-type m))
         (zero (coerce 0 element-type))
         (pos-k (coerce (/ alpha lambda) element-type))
         (neg-k (- pos-k)))
    (dotimes (i (array-total-size m))
      (let* ((m-i (row-major-aref m i))
             (v-i (cond ((> m-i pos-k)
                         (- m-i pos-k))
                        ((< m-i neg-k)
                         (+ m-i neg-k))
                        (t
                         zero))))
        (setf (row-major-aref m i) v-i))))
  (values))

(defun nsolve/double-float (alpha lambda m)
  (declare (type double-float alpha lambda)
           (type (simple-array double-float (*)) m)
           (optimize speed))
  (%nsolve alpha lambda m))

(defun nsolve/default (alpha lambda m)
  (declare (type real alpha lambda)
           (type array m))
  (%nsolve alpha lambda m))

The runtime performance of the nsolve function which is specialized to double floats is

;; Timing: nsolve/double-float
;; Evaluation took:
;;   0.224 seconds of real time
;;   0.223517 seconds of total run time (0.222565 user, 0.000952 system)
;;   100.00% CPU
;;   448,054,948 processor cycles
;;   0 bytes consed

The reduction in time is obviously preferable. Also note that the amount of consing has reduced to zero.

Unfortunately, improving performance has introduced a number of other problems. We have lost the generic interface to the nsolve function as the user must select the most appropriate specialization manually.

The implementor of the proximal operator could solve this problem by providing an nsolve function which dispatches to the most applicable specialization. This would burden the developer with the maintenance of the specializations, the dispatch function and ensuring portability across common lisp implementations.

The maintenance burden would be reasonable if there were only one operator, but there are many proximal operators

• l2 proximal operator
• l2 squared proximal operator
• l-infinity proximal operator
• Nuclear norm proximal operator
• Spectral norm proximal operator
• inequality proximal operator
• subspace proximal operator
• ...

Furthermore, these operators are also used inside optimization strategies like Accelerated Proximal Gradient or the Alternating Direction Method of Multipliers.

All of these problems would require a function that can be inlined (a lambda form), a set of specializations and a dispatch function (runtime and/or compile time).

The purpose of the template function system is to simplify the implementation of such problems. The system is responsible for inserting specializations in to the global environment and the management of the dispatch function.

The template function system does not create specializations on demand. A user must request a specialization prior to its use. This strategy shifts the burden of portability across common lisp implementations to the application. This makes sense since it is the application which imposes the requirements (restrictions) on representation, not the operator.

We hope that this example sufficiently motivates the need for the template function system. The system should find many uses in operators that perform mathematical optimization and operators that are defined for sequences.