# Generating an ordered list of pairs of elements from ordered lists

I have a pair of ordered lists. I want to generate a new ordered list (using the same ordering) of length n by applying a binary operator to pairs of elements, one from each list, along with the index of each list's elements.

For example, if I have two lists of real numbers in increasing order (sorted via Less),

l1 = {5.7832, 30.4713, 74.887, 139.04, 222.932, 326.563}

l2 = {3.481, 9.2816, 15.7112, 27.8226, 45.0379, 67.222, 94.3526}


and the binary operator is Plus, then if I ask for 5 elements, the result would be

{{9.26417, 1, 1}, {15.0647, 1, 2}, {21.4944, 1, 3}, {33.6058, 1, 4}, {33.9522, 2, 1}}


I've written code that achieves this but I'm looking for the most elegant, efficient, and general solution that would work with arbitrary lists and binary operators.

-
By the way, welcome to Mathematica.SE! For some reason I thought you had already joined, early on, but can't find the account now in the user list. –  Verbeia Sep 11 '12 at 10:34
Hey Paul, nice to see you here. We met last year at the WRI tech conference (though you probably won't remember that). –  Sjoerd C. de Vries Sep 11 '12 at 10:40
I notice that all responders so far have solved the problem, but I wonder about the efficiency aspect. The solutions generally generate all tuples, use the operator and then sort, and pick the first few elements. This is fine for short lists, but with n=10^6 it will be difficult, whereas it should be possible to pick out the first 5 elements without bruteforcing as it is given that the lists are sorted. –  Sjoerd C. de Vries Sep 11 '12 at 11:38
@SjoerdC.deVries For a generic binary operator (i.e. non monotonic) how to use the information about the input lists being sorted ? –  b.gatessucks Sep 11 '12 at 13:00
Ahoi Paul, welcome to the party! –  Yves Klett Sep 11 '12 at 13:42

As a variation of b.gatessucks's answer:

orderedPairs[oper_, list1_List, list2_List, n_Integer?Positive] :=
Take[Sort[Flatten[MapIndexed[Flatten[{##}] &, Outer[oper, list1, list2], {2}], 1]], n]

-

## Preamble

I can't resist to answer, although the answer won't be short :-). Below is an attempt to systematically capitalize on properties which weren't used (or at least used systematically) in other answers.

I will use the same lists provided by the OP, for tests:

l1 = {5.7832, 30.4713, 74.887, 139.04, 222.932, 326.563}

l2 = {3.481, 9.2816, 15.7112, 27.8226, 45.0379, 67.222, 94.3526}


## The algorithm

The main idea is that for many common binary operators, the order of the list resulting from applying the operator to a list and a single element of another list can be deduced from the order of the original list and the element in question, often without explicit sorting of the resulting list.

The algorithm itself will be based on lazy merging of lazy lists, rsulting from application of the binary operator in question to a list (say l1) and each of the elements of the other list l2. What we will do is to create a binary search tree,which holds the information nodes of the following form:

{{currentNumber, index1,index2},iterator}


where an iterator is a construct which supplies the next element of a resulting list on demand. What is important is that an iterator encodes the information about the ordering of the new sub-list.

What we will do is to merge the resulting sub-lists lazily, very analogously to what merge function does for 2 lists in the mergesort algorithm, but now we have Length[l2] sub-lists (represented by their iterators). Here is the outline of the general algorithm to get a single element with indices:

1. Create Length[l2] iterators, each encoding a lazy version of the result of an application of our binary operator to a list l1 and i-th element of l2.

2. Construct an initial binary search tree by requesting each iterator to provide the corresponding first element and constructing nodes of the form outlined above. The ordering used for the binary seacrh tree construction can be also specified.

3. If we already picked enough elements, exit.

4. Pick the node with the smallest (according to our ordering) element, and save it in a temporary variable.

5. Delete this node from a tree.

6. Construct the new node by requesting the next element from the iterator of the smallest (deleted) node, if that iterator still has elements, and add this node to the tree.

7. Return the first part of the deleted node, {currentNumber, index1,index2}

8. Go to 3

Since insertion and deletion time for a (balanced) binary search tree is O(Log(m)), where m is the number of nodes (which is the current number of iterators, can be at most Length[l2]), the complexity of the search alone will be of the order of n Log[m], where n is the number of elements requested.

This can be very fast for small enough n and large enough lists, but the main question is whether or not we can build the iterators efficiently. In general, this would involve copying and sorting all sub-lists, and therefore should have a complexity Length[l1]Log[Length[l1]]Length[l2]. However, for some common cases (such as Times andPlus), the resulting ordering can be deduced from the initial ordering and the element of the second list in just constant time, and then we win big both speed-wise and memory-wise.

## Mathematica impelentation

### Binary search tree

For simplicity, I will use the unbalanced binary search tree implementation, which is a slightly generalized (to allow a custom ordering function) version of an excellent implementation due to Roman Maeder ("Computer Science With Mathematica", chapter 6). Here it is:

These are constructors and selectors for the tree nodes:

Clear[node, tree, leftTree, rightTree, emptyTree];
node[info_, left_, right_] := tree[info, left, right]
leftTree[tree[_, left_, _]] := left
rightTree[tree[_, _, right_]] := right
information[tree[info_, _, _]] := info


This is the insertion function:

Clear[insertTree];
insertTree[emptyTree, info_, Order_: Order, Key_: Identity] :=
node[info, emptyTree, emptyTree]

insertTree[tree_, info_, Order_: Order, Key_: Identity] /;
Order[Key[info], Key[information[tree]]] == 0 := tree;

insertTree[tree_, info_, Order_: Order, Key_: Identity] /;
Order[Key[info], Key[information[tree]]] > 0 :=
node[
information[tree],
insertTree[leftTree[tree], info, Order, Key],
rightTree[tree]
];

insertTree[tree_, info_, Order_: Order, Key_: Identity] :=
node[
information[tree],
leftTree[tree],
insertTree[rightTree[tree], info, Order, Key]
]


Node insertion is a relatively simple matter since one just needs to find the place where to insert it, and the way the tree is constructed automatically preserves the right order.

We won't need the node search functionality, so I skip it. This is the function to find the smallest node by straightforward tree descent:

Clear[smallestNode];
smallestNode[emptyTree] = Null
smallestNode[tree_] /; leftTree[tree] === emptyTree :=
information[tree]
smallestNode[tree_] := smallestNode[leftTree[tree]]


Here is the most conceptually difficult function to delete the node:

Clear[deleteTree];
deleteTree[tree_, key_, Order_: Order, Key_: Identity] /;
Order[Key[information[tree]], key] < 0 :=
node[
information[tree],
deleteTree[leftTree[tree], key, Order, Key],
rightTree[tree]
]

deleteTree[tree_, key_, Order_: Order, Key_: Identity] /;
Order[Key[information[tree]], key] > 0 :=
node[
information[tree],
leftTree[tree],
deleteTree[rightTree[tree], key, Order, Key]
]

deleteTree[tree_, key_, Order_: Order, Key_: Identity] /;
leftTree[tree] == emptyTree := rightTree[tree]

deleteTree[tree_, key_, Order_: Order, Key_: Identity] /;
rightTree[tree] == emptyTree := leftTree[tree]

deleteTree[tree_, key_, Order_: Order, Key_: Identity] :=
With[{nextinfo = smallestNode[rightTree[tree]]},
node[
nextinfo,
leftTree[tree],
deleteTree[rightTree[tree], Key[nextinfo], Order, Key]
]
]


Finally, here is the constructor to construct tree from a list of information items:

Clear[Tree];
Tree[l_List, Order_: Order, Key_: Identity] :=
Fold[Function[{t, i}, insertTree[t, i, Order, Key]], emptyTree, l]


Some examples:

(t = Tree[{5, 4, 7, 3, 2, 8, 6}]) // TreeForm
deleteTree[t, 5] // TreeForm


I may include some visualization function to plot a tree better than TreeForm, in the future.

### Iterators

Here is a generic iterator "data type":

Clear[getNext, iterator];
getNext[iterator[fun_]] := fun[];
iterator /: MakeBoxes[it_iterator, fmt_] :=
InterpretationBox[RowBox[{"«", "iterator", "»"}], it];


It basically tells that objects with heads iterator contain a function embedded in them, and the end user should only know that applying to the getNext, s/he gets the next element.

We now become specific and construct the iterator generator for iterators for Plus binary operator:

Clear[makePlusInterator];
makePlusInterator[l1_, elem_] :=
With[{len = Length[l1]},
Module[{i = 0},
iterator[Function[If[i < len, l1[[++i]] + elem, None]]]
]
];

makePlusInterator[l1_, l2_List] :=
Map[makePlusInterator[l1, #] &, l2];


Our iterators so constructed are closures, embedding addition as a part of producing an element on demand. Here we encode our knowledge of the resulting ordering (which, in the case of Plus, is the same as the original one).

Here is an example:

iter = makePlusInterator[l2,First@l1]

(*  « iterator »  *)

Table[getNext[iter],{6}]

(*  {9.2642,15.0648,21.4944,33.6058,50.8211,73.0052}  *)


Note that iterators themselves do not carry the index information (perhaps it would be more logical if they would).

### Bringing it together

To combine the components I just described and actually implement the algorithm, we need two functions. First is a helper function which implements a part of the step 6:

Clear[getNextInfo];
getNextInfo[{{next_, fstind_, secind_}, iter_}] :=
{{getNext[iter], fstind, secind + 1}, iter};


It takes the information stored in a node and prepares the information for the new node to be inserted into the tree. Here is the main function:

Clear[makeResultinTreeInterator];
makeResultinTreeInterator[iterators : {__iterator}, Order_: Order] :=
Module[{
tr =
Tree[
MapIndexed[{{getNext[#], First@#2, 1}, #} &, iterators],
Order,
First
]
},
iterator[
Function[
Module[{smallest = smallestNode[tr], next},
tr = deleteTree[tr, smallest];
next = getNextInfo[smallest];
If[First@next =!= None,
tr = insertTree[tr, next]
];
First@smallest]
]
]
];


Note that the function itself returns an iterator, so our resulting list is lazy as well.

### Examples

I will now illustrate this:

resiter = makeResultinTreeInterator[makePlusInterator[l2, l1], Order]
Table[getNext[resiter],{6}]
Table[getNext[resiter],{6}]
Table[getNext[resiter],{6}]

(*

{{9.2642,1,1},{15.0648,1,2},{21.4944,1,3},{33.6058,1,4},{33.9523,2,1},{39.7529,2,2}}

{{46.1825,2,3},{50.8211,1,5},{58.2939,2,4},{73.0052,1,6},{75.5092,2,5},{78.368,3,1}}

{{84.1686,3,2},{90.5982,3,3},{97.6933,2,6},{100.136,1,7},{102.71,3,4},{119.925,3,5}}

*)


So far so good. Now, we pick some larger lists:

large1 = Sort@RandomReal[{-100, 100}, 1000];
large2  = Sort@RandomReal[{-100, 100}, 100];
randomizedOrderLarge2 = RandomSample[large2];
resiterLarge =
makeResultinTreeInterator[
makePlusInterator[large1 , randomizedOrderLarge2], Order];


Here I had to randomize the second sample because the implementation of the binary search tree I use is for unbalanced trees, which would degenerate to a line for an ordered second list. This somewhat violates the initial condition that both lists are ordered, but I only use this for benchmarks, and the resulting list won't change in any case.

We now test by fetching some larger number of elements:

Table[getNext[resiterLarge],{500}]//Short//AbsoluteTiming
Table[getNext[resiterLarge],{1000}]//Short//AbsoluteTiming
Table[getNext[resiterLarge],{2000}]//Short//AbsoluteTiming

(*
{0.3925781,{{-191.248,33,1},<<498>>,{-172.676,27,43}}}
{0.8896485,{{-172.665,6,53},<<998>>,{-161.276,29,162}}}
{1.9609375,{{-161.251,56,14},<<1998>>,{-144.795,20,134}}}
*)


We see that the timing linearly depends on the number of elements, but our code is really slow.

## Java implementation

Now that we have our slow Mathematica prototype, we can port it to Java in the hopes that this will run (much) faster. To work with this part, you will need the Java reloader (which I modified to recognize interfaces, so if you have the old one you should change to use the recent version).

While what follows is conceptually a direct port of the Mathematica part, there will be significant differences in details. In particular, I will use the Java TreeSet library class, so that I don't need to implement the binary search tree from scratch.

So, run the Java reloader code first.

### General interface for customization

The following interface should be implemented for each binary operator and ordering, to customize the algorithm (execute this code):

JCompileLoad@
"package orderedPairsFinal;

import java.util.Comparator;
import java.util.Iterator;

public interface Customizer {
public Iterator<Double> getIterator(double[] list, double elem);
public Comparator<Double> getComparator();
}"


It basically declares that we have a way to compare doubles and construct iterators (I stick to doubles here).

### Customization for Plus

This implements the above interface for Plus:

JCompileLoad@"
package orderedPairsFinal;

import java.util.Comparator;
import java.util.Iterator;

public class  PlusCustomizer implements Customizer {
@Override
public Comparator<Double> getComparator() {
return new Comparator<Double>() {
public int compare(Double arg0, Double arg1) {
return Double.compare(arg0, arg1);
}
};
}

@Override
public Iterator<Double> getIterator(final double[] list,
final double elem) {
return new Iterator<Double>() {
private int ctr = 0;

public boolean hasNext() {
return ctr < list.length - 1;
}

public Double next() {
return list[ctr++] + elem;
}

public void remove() { }

};
}
}"


### Tree node class

The following class implements the tree node abstraction.

JCompileLoad@"
package orderedPairsFinal;

import java.util.Comparator;
import java.util.Iterator;

public class NodeInfo implements Iterator<NodeInfo>,Comparable<NodeInfo> {
private Iterator<Double> iter;
private Comparator<NodeInfo> comp;
public Double current;
public final int fstIndex;
public int secIndex;

public NodeInfo(Iterator<Double> iter,Comparator<NodeInfo> comp,
int index, int secIndex) {
this.fstIndex = index;
this.secIndex = secIndex;
this.iter = iter;
this.comp = comp;
current = iter.next();//Note: no checks performed here
}

@Override
public NodeInfo next() {
return new NodeInfo(iter,comp,fstIndex,secIndex+1);
}

@Override
public boolean hasNext() {
return iter.hasNext();
}

public void remove() { }

@Override
public int compareTo(NodeInfo arg0) {
return comp.compare(this, arg0);
}

}"


### The main class

This is the class which implements our algorithm.

JCompileLoad@"package orderedPairsFinal;

import java.util.Comparator;
import java.util.TreeSet;

public class OrderedPairs {
private double[] first;
private double[] second;
private Comparator<NodeInfo> comp;
private TreeSet<NodeInfo> tree = new TreeSet<NodeInfo>(comp);

public OrderedPairs(Customizer cust, double[] list1,
double[] list2) {
first = list1.length>list2.length?list1:list2;
second = list1.length>list2.length?list2:list1;
final Comparator<Double> dcomp = cust.getComparator();
comp =
new Comparator<NodeInfo>() {
public int compare(NodeInfo o1, NodeInfo o2) {
return dcomp.compare(o1.current, o2.current);
}
};
for(int i=0;i<second.length;i++){
}
}

public double[] getNext(){
if(tree.isEmpty()){
return null;
}
NodeInfo fst = tree.first();
tree.remove(fst);
if(fst.hasNext()){
}
return new double[]{fst.current, fst.fstIndex+1,fst.secIndex+1};
}

public double[][] getNext(int num){
double[][] aux = new double[num][3];
int ctr = 0;
while(ctr < num){
double[] next = getNext();
if(next == null) break;
aux[ctr++] = next;
}
double[][] res = new double[ctr][3];
System.arraycopy(aux, 0, res, 0, ctr);
return res;
}
} "


The methods of this class and the way it works require some explanation. The class takes a customizer and two arrays of doubles in constructor, and creates iterators and the tree also in constructor. It has the method getNext, which is overloaded to have 2 forms. You can either request a single element (with indices), or do that in a bulk. The latter possibility is particularly relevant for uses from Mathematica, to avoid the killing overhead of many individual method calls in JLink.

### Examples

To use this, call:

cust = JavaNew["orderedPairsFinal.PlusCustomizer"];
op  = JavaNew["orderedPairsFinal.OrderedPairs", cust, l1, l2];


We can now use it:

op@getNext[6]
op@getNext[6]

(*
{{9.2642,1.,1.},{15.0648,1.,2.},{21.4944,1.,3.},{33.6058,1.,4.},
{33.9523,2.,1.},{39.7529,2.,2.}}

{{46.1825,2.,3.},{50.8211,1.,5.},{58.2939,2.,4.},{73.0052,1.,6.},
{75.5092,2.,5.},{78.368,3.,1.}}
*)


Now, the larger benchmarks:

large1 = Sort@RandomReal[{-100, 100}, 1000];
large2  = Sort@RandomReal[{-100, 100}, 1000];
opL  = JavaNew["orderedPairsFinal.OrderedPairs", cust, large1, large2];


and we can test:

opL@getNext[50000]//Short//AbsoluteTiming
opL@getNext[50000]//Short//AbsoluteTiming

(*
{0.0732421,{{-199.219,1.,1.},<<49998>>,{-137.102,73.,243.}}}
{0.0419922,{{-137.102,17.,304.},<<49998>>,{-110.502,335.,108.}}}
*)


which is quite fast for lazy lists to my mind, and this time includes the data transfer!(possibly being dominated by it).

## Conclusions

I presented the abstract algorithm which promises good overall performance and very good performance for practically important special cases. I have also presented two implementations of it, using Mathematica and Java.

While Mathematica implementation has the right algorithmic behavior, it is very slow. This is one of the cases where the symbolic overhead of full Mathematica evaluator really gets in the way. The data structures involved and the structuring of code which allows to separate customizable parts from general parts (iterators, lazy lists, binary search tree) here did not allow to easily Compile any part of the algorithm, so there are little chances to improve its speed within Mathematica.

Implementation in Java allowed to use powerful built-in data structures and idioms which in Java do not cost such an overhead as in Mathematica. As a result, the Java version is thousands times faster.

In both cases, I could also include the customizations for other common binary operators such as Times, but did not do so simply because this post is already excessively long.

-
+1 Definite food for thought re Java OO integration. V.general implementation. Would it be difficult to integrate optimisations that capitalise on the structure of the two lists? E.g. by looking ahead in l2 and noting that p=27.8226 is the largest element of this list such that 5.7832+p <3.481+30.4713 and therefore getting all pairs preceding {33.6058, 1, 4} “for free”. l1 = {5.7832, 30.4713, 74.887, 139.04, 222.932, 326.563} l2 = {3.481, 9.2816, 15.7112, 27.8226, 45.0379, 67.222, 94.3526} {{9.26417, 1, 1}, {15.0647, 1, 2}, {21.4944, 1, 3}, {33.6058, 1, 4}, {33.9522, 2, 1}} –  Ronald Monson Sep 12 '12 at 22:47
@RonaldMonson Thanks for the upvote. As to your suggestion, I will think about it, but my first impression is that it is hard to do such optimizations generally. One can probably construct more complex scheme of merging, where when we know that we should take further values from the same node / iterator and if the node/ interator remains in the same place in the tree after we extract one or more elements, we can skip repeated tree searches / insertions and take values from the iterator directly. But even now, the algorithm is already quite fast, and it is lazy as well. So, basically, you ... –  Leonid Shifrin Sep 12 '12 at 22:52
@RonaldMonson ... get a single extra element on demand in quazi-constant (multiplied by a log of the length of one of the smallest of the two lists) time for many operators, what can be much better than that? Of course, some time is spent because I use some additional objects to make things flexible. I can imagine a C version with ahead-of-time compilation similarly in spirit to what R.M does in his answer, with yet much better (than the current one of my answer) performance. –  Leonid Shifrin Sep 12 '12 at 22:54
Yes, it perhaps unlikely to be materially different in this particular case but I'm wondering in the general; Mma is great for testing optimizations but is it a case of diminishing efficiency gains as it becomes increasingly difficult to implement refinements in the corresponding Java version? –  Ronald Monson Sep 12 '12 at 23:06
@RonaldMonson I would think that in a case like this, there is a high chance that efficiency gains can be hidden in the mma version by a huge symbolic overhead. OTOH, since my mma implementation has (hopefully) correct algorithmic behavior, and the overhead is just in the huge constant multiplicative factor, one should be able to comparatively analyze modifications in mma. As to Java, the current algorithm is quite simple (you can see that the algorithm itself is only a dozen lines or so), and it should not be hard to modify that one as well. Java development is harder for a different ... –  Leonid Shifrin Sep 12 '12 at 23:14

The solutions using Outer and Tuples will all blow up for large lists. The following is a compiled function that exploits the ordered nature of the lists and should be very fast and low on memory footprint even on very large lists.

cf = With[{op = Plus},
Compile[{{list1, _Real, 1}, {list2, _Real, 1}, {n, _Integer}},
Module[{result = InternalBag[], count = 0, i, j,
len1 = Length@list1, len2 = Length@list2, val, uBound},

For[i = 1, i < len1, i++,
uBound = Max[op[list1[[i + 1]], list2[[i]]], op[list1[[i]], list2[[i + 1]]]];
If[count < n,
For[j = 1, j < len2, j++,
If[(val = op[list1[[i]], list2[[j]]]) <= uBound && count < n,
count += 1;
InternalStuffBag[result, {val, i, j}, 1];,
Break[]
]
],
Break[]
]
];
InternalBagPart[result, All] ~Partition~ 3
],
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
]
]


You will have to create a different compiled function for a different binary operator. You can probably even write a wrapper around this to generate compiled code on the fly (assuming your binary operator is compilable).

Using the above function on your lists:

cf[l1, l2, 5]
(* {{9.2642, 1., 1.}, {15.0648, 1., 2.}, {21.4944, 1., 3.}, {33.6058, 1., 4.}, {33.9523, 2., 1.}} *)

-
+1 This is what I meant with my comment above. I'll pretend I didn't see the For-s ;-) –  Sjoerd C. de Vries Sep 11 '12 at 19:04
Any reason why you chose For[] instead of Do[]? :) –  Guess who it is. Sep 12 '12 at 2:04
@J.M. Lol, no particular reason other than it doesn't matter in compiled code and this might be simpler to follow for someone reading procedural code :) Boy, I've gotten a few comments about this already here and in chat... ;) –  The Toad Sep 12 '12 at 2:12
Gives errors in the following case: l1 = Sort@RandomReal[{-100, 100}, 100];l2 = Sort@RandomReal[{-100, 100}, 10]; cf[large1, large2, 100]. Generally, your code seems to assume the monotonous behavior of the total resulting ordering, which is true for Plus, but not true in general, e.g. for Times when lists contain negative numbers (or may be I misunderstood). Very nice code still, I voted for it shortly after you posted. –  Leonid Shifrin Sep 12 '12 at 15:01
@LeonidShifrin Ah, in that case you must have used the smaller list as the first argument... This happens when it tries to index the smaller list with a number larger than its length (happens when n is comparable to the length or larger). I'm aware of this and I should've added a few additional checks and a more robust way of getting the bound and doing the ordering. I thought it would be as simple as flipping the order of the lists, but I now see it isn't the case (at least, not yet). I'll get back to this later today and fix some of the issues –  The Toad Sep 12 '12 at 16:32

There is another way using RankedMin (we don't have to sort the whole list as we would do with Sort, SortBy) :

RankedOrd[l1_, l2_, n_, f_] :=
With[{l = Flatten[Array[{f[l1[[#1]], l2[[#2]]], #1, #2} &, Length /@ {l1, l2}], 1]},
l[[#]] & /@ Table[Position[l, RankedMin[Transpose[l][[1]], i]][[1, 1]], {i, n}]]


for the lists l1, l2 and the function Plus we have :

RankedOrd[l1, l2, 5, Plus]

{{9.2642, 1, 1}, {15.0648, 1, 2}, {21.4944, 1, 3}, {33.6058, 1, 4}, {33.9523, 2, 1}}

-
Did you use the same l1 and l2 as in the example? The output differs completely. –  Sjoerd C. de Vries Sep 11 '12 at 11:25
Yes, but I assumed he wanted 5 max elements instead of 5 min. In that case I'd use [[;;5]] instead of // Take[#,-5] & as well as RankedMin. –  Artes Sep 11 '12 at 11:30
@SjoerdC.deVries I changed the solutions to get the output as the OP wanted. –  Artes Sep 11 '12 at 15:39

Something like :

f[oper_, list1_, list2_, n_] := Sort[Flatten[
Outer[{oper[list1[[#1]], list2[[#2]]], #1, #2} &,
Range[Length[list1]], Range[Length[list2]]], 1]][[1 ;; n]]

f[Plus, l1, l2, 5]

(* {{9.2642, 1, 1}, {15.0648, 1, 2}, {21.4944, 1, 3}, {33.6058, 1,  4}, {33.9523, 2, 1}} *)


Also :

f2[oper_, list1_, list2_, n_] := Sort[{oper[#[[1, 1]], #[[2, 1]]], #[[1, 2]], #[[2, 2]]} & /@
Tuples[{Transpose[{list1, Range[Length[list1]]}],
Transpose[{list2, Range[Length[list2]]}]}]][[1 ;; n]]

f2[Plus, l1, l2, 5]

(* {{9.2642, 1, 1}, {15.0648, 1, 2}, {21.4944, 1, 3}, {33.6058, 1, 4}, {33.9523, 2, 1}} *)


In the specific case that l1, l2 are sorted and the operator is monotonic then the wanted first n results must come from using the first n elements of the input lists. In this case the functions can be used more efficiently as :

f[Plus, l1[[1 ;; 5]], l2[[1 ;; 5]], 5]
f2[Plus, l1[[1 ;; 5]], l2[[1 ;; 5]], 5]

-
Interestingly, the 2nd and last number in the output list have a different last decimal than the one in the given example. Apparently MachinePrecision is rather different on both computers or so? –  Sjoerd C. de Vries Sep 11 '12 at 10:48
@SjoerdC.deVries Wow, sharp eye ! I don't know, it would be useful to know how the output was produced. –  b.gatessucks Sep 11 '12 at 10:52

Another variation:

Sort[{Array[{l1[[#1]] + l2[[#2]], ##} &, Length /@ {l1, l2}, {1, 1}, Sequence]}][[;; 5]]

-

This answer contains no Mma code.

Naturally there will not be a single solution that is the most “elegant, efficient or general” so either one could produce 3 champions on each of these categories or a single solution that has the right trade-offs (in which case perhaps further criteria would be needed) - personally I’d also include “comprehensibility” since I think this is sufficiently distinct from “elegant” to be useful.

None of the provided solutions to date offer a completely general solution since although arbitrary lists and binary operators are specified in the question, the key factor in the complexity and implementation of this problem is the Sort function; in particular, the relationship between the OrderedQ functions used to sort both initial lists (OrderedQ1,OrderedQ2) and the output of the binary operator (OrderedQ3).

If these OrderedQ functions obey a type of transitivity in relation to the operator op (essentially we can use the original sorting in the final sort):

(T1)

1. $\forall$ a in l1, $\forall$ b,c in l2, OrderedQ2[b,c] => OrderedQ3[op[a,b],op[a,c]]
2. $\forall$ x in l2, $\forall$ y,z in l1, OrderedQ1[y,z] => OrderedQ3[op[x,y],op[x,z]]

then only the first n elements of both lists need be considered with the resulting function having complexity $O(n^2)$. Clearly in this example T1 is satisfied with Plus and Less so that n up to a thousand should be manageable (with no restriction on list size since elements beyond the nth position can be thrown away)

On the other hand if T1 is not satisfied then a worse-case scenario would involve comparing all pairs from both lists (of size say $m$) before sorting the result and taking the first n elements- a function with complexity $O(m^4)$ and only lists in the hundreds/low thousands could be managed.

In this non-T1 scenario, improvements in the worst-case complexity might come from considering the relationship between the three OrderedQ functions, and the structure/ distribution within and between l1 and l2. Any regularities could probably be tackled in a general way with by recursing down both l1 and l2 and/or using using the problem’s inherent parallellizability.

In summary:

• With T1: R.M’s compiled function will be quickest (but not as readable) and all the other answers would get galactic speed-up by migrating the "1;;5" specification inside Sort thereby only considering the first n elements (Artes's solution only sorts the first n of the binary-operator output by using RankedMin but still applies the binary operator to all pairs) .
• Without T1: R.M’s compiled function is no longer applicable (even if the binary operator is compilable) and all the other solutions are easily modifiable by adding the appropriate OrderedQ3 to Sort but with a significant performance hit (in comparison with the aforementioned migration).

• More subjectively; belisarius’s solution is the most elegant IMO.

-
I don't really follow what you're trying to say here... The problem with the other answers is not in where the 1;;5 portion lies, but rather that Outer and Tuples will generate gigantic lists that'll quickly blow up memory. If you're saying take only the first 5, then you ought to consider that 5 here is just an example... it might as well be 10^5. Secondly, the compiled solution can be easily modified to use a different ordering although, it must be recompiled and the ordering function must be written such that it is compilable –  The Toad Sep 12 '12 at 6:46
Granted, the compiled solution comes with a small loss in readability (especially if you're alien to reading procedural code) and for reasonably sized lists, I would probably use Outer myself –  The Toad Sep 12 '12 at 6:48
My point was that most of the work in Outer and Tuples is not needed if the size of the two lists is >>n since you only need to worry about the first n elements ( a fact that actually only the compiled version seems to have taken into account - provided T1 is accepted). Without T1 however the whole lists needs to be accounted for unlike the compiled version (if I understand its algorithm correctly -yes it is a little alien to me). –  Ronald Monson Sep 12 '12 at 6:51
Could you explain what you mean by your T1 – maybe I didn't interpret it properly... In the problem presented, OrderedQ1 = OrderedQ2 = OrderedQ3, i.e. the same ordering is used for both lists and n elements of the output are desired following the same ordering. I took his "arbitrary list" to mean any list (albeit sorted) of unknown length, whereas I think you take it to mean arbitrary sorting. Is that a correct interpretation? –  The Toad Sep 12 '12 at 6:57
Yes. I actually edited out the following in the original answer as I may be being a little pedantic here - "(OK helps to read the question - I think T1 is assumed from "using the same ordering" but then the binary operator is implicitly restricted to have same domain range and all solutions except - R.M' s have much higher complexity than necessary by considering elements beyond the nth position) (No, "using the same ordering" doesn' t necessarily imply T1 so the question was read properly :))" - Is this right, does your compiled version only use the first n values and ignores the rest? –  Ronald Monson Sep 12 '12 at 7:03
ClearAll[ord, ordPairs];
ord[list_List, comp_, oper_, n_] :=
With[{lst = oper /@ list}, list[[Ordering[lst, n, comp]]]];
ordPairs[list1_List, list2_List, comp_, oper_, n_] :=
With[{ls1 = ord[list1, comp, Identity, All],
ls2 = ord[list2, comp, Identity, All]},
{oper[#], First @@ Position[ls1, #[[1]]], First @@ Position[ls2, #[[2]]]} & /@
ord[Tuples[{ls1, ls2}], comp, oper, n]]


Examples:

ordPairs[l1, l2, Less, Total, 5]
(*{{9.2642, 1, 1}, {15.0648, 1, 2}, {21.4944, 1, 3}, {33.6058, 1, 4}, {33.9523, 2, 1}}  *)
ordPairs[l1, l2, Greater, Total, 5]
(* {{420.916, 1, 1}, {393.785, 1, 2}, {371.601, 1, 3}, {354.386, 1, 4}, {342.274, 1, 5}}*)
ordPairs[l1, l2, Greater, Times @@ # &, 3]
(* {{30812.1, 1, 1}, {21952.2, 1, 2}, {21034.2, 2, 1}}*)
ordPairs[l1, l2, Less, Times @@ # &, 3]
(* {{20.1313, 1, 1}, {53.6773, 1, 2}, {90.861, 1, 3}}*)
ordPairs[l1, l2, Less, Plus @@ # &, {3}]
(* {{21.4944, 1, 3}} *)
ordPairs[l1, l2, Less, Plus @@ # &, {-3}]
(* {{371.601, 6, 5}} *)

-

Golfing:

Sort@Transpose[{Total /@ Tuples@#, Tuples @@ {Range /@ Length /@ #}}][[;; 5]] &@{l1, l2}

-

To clarify the comment I posted below the question:

If I have an binary operation f with property

$f(a_1,b_1) \geq f(a_2,b_2)$ if $a_2\geq a_1 \wedge b_2\geq b_1$

(examples of such operators are Times and Plus), and if I have a sorted list $a_1,..,a_{n+1}$ and $b_1,..,b_{n+1}$ then it must hold that

$f(a_{n+1},b_1)\geq f(a_{i},b_1)\wedge f(a_1,b_{n+1})\geq f(a_1,b_i)$ for any $1\leq i\leq n$,

because due to the lists being sorted $a_{n+1}\geq a_i$ and $b_{n+1}\geq b_i$. Hence, there must be at least $n$ results ($f(a_1,b_1), f(a_2,b_1),..,f(a_n,b_1)$ or $f(a_1,b_1), f(a_1,b_2),..,f(a_1,b_n)$) that are smaller than what can obtained using either of the $n+1$th elements.

I realize that this explanation overlaps with that of Ronald, but he sort of sidesteps the proof, by simply asserting that "then only the first n elements of both lists need be considered".

-
MapIndexed[#2~Prepend~#1 &,
Outer[Plus, l1, l2], {2}]~Flatten~1~SortBy~First


or very similarly

FlattenAt[
MapIndexed[List, Outer[Plus, l1, l2], {2}]~Flatten~1~SortBy~First,
Tuples@{Range[Length@l1 Length@l2], {2}}]


or

DeleteCases[
MapIndexed[List, Outer[Plus, l1, l2], {2}]~Flatten~1~SortBy~
First, List, {3}, Heads -> True]


or

Join[Map[List, Outer[Plus, l1, l2], {2}],
Outer[List, Range@Length@l1, Range@Length@l2],
3]~Flatten~1~SortBy~First

-