Formulation of the Problem

What are the goals of science and engineering. We want to know what will happen in the future, and we want to decide what to do to make this future better. In a nutshell:

• predicting the future is what science is about, and

• deciding what to do is more the area of engineering.

Need for measurements and computations. To predict the world's future and to decide how to improve this future, we need to understand how the world works. Some of the situations we can simply observe, but in most cases, we need to measure the values of different quantities.

Based on the measurement results, we deduce the algorithms describing dynamics of the world's processes -and formulate selection of the best action and/or device as an appropriate optimization problem. To use prediction algorithms and to solve the corresponding optimization problem, we often need to perform a large number of computations. For example, to reasonably accurately predict tomorrow's weather, we need to use high-performance computers.

Comment. Some problems require so many computational steps that we cannot solve them even on the fastest computers. For example, to predict in which direction a tornado will go, even the existing high-performance computers are not sufficient:

ICT&ES-2020: Information-Communication Technologies & Embedded Systems, November 12, 2020, Mykolaiv, Ukraine vladik@utep.edu (V. Kreinovich); mceberio@utep.edu (M. Ceberio); olgak@utep.edu (O. Kosheleva) 0000-0002-1244-1650 (V. Kreinovich); 0000-0003-2587-4209 (O. Kosheleva) • to predict where the tornado will do in the next 15 minutes, we need to spend several hours on a high-performance computer;

• by that time, the tornado has already moved and caused damage.

Our resources are limited. In principle, for the same computational problem, we may have different algorithms. Which algorithm should we select? Traditional approach of theoretical computer science uses asymptotic analysis to compare different algorithms. For example, we analyze how the computation time grows with the size 𝑛 of the input, and we select an algorithm that leads to the shortest possible time for large 𝑛. This analysis implicitly assumes that we can have an unlimited amount of resources:

• we can have inputs of arbitrary large size 𝑛,

• we can have an arbitrary large memory to store this data and to store all needed intermediate results,

• we can perform an arbitrary large number of computations, etc.

In practice, our resources are limited. It is therefore important to consider how this limitation affect computations and measurements.

What we do in this paper. In this paper:

• we provide a brief overview of the corresponding resource limitations, and

• we also present two case studies explaining, in some detail, how to take these resource limitations into account.

Comment. In most cases, we focus on our own related results -although, of course, in each of these problems, there are many other results and techniques. The reason for this selection is that the results that we cite are usually mathematically reasonably simple -this is why we selected them: for our own results, we know how we can simplify these results into easy-to-follow not-too-complex ones. We believe that overall, these results provide a good introduction to the great variety of related problems.

White-box and black-box computing.

In our analysis, we will distinguish between:

• traditional (white-box) computing and • practically important black-box computing.

White-box computing corresponds to the usual analysis of theoretical computer science, when we know, step-by-step, what our algorithms are doing. This corresponds to:

• situations when we write our own code "from scratch", and

• situations when we only use open-source software.

In practice, however, we often use proprietary code, i.e., a code that its producer does not explain. We can use it, but we are not allowed to see the algorithm -we only see the results. Similarly, in defense applications, we may need to use classified code -e.g., the code used to control the corresponding devices.

White-Box Computing: Three Types of Situations and Related Resource Limitations

Three types of situations. In our analysis, we will distinguish between three different computational situations.

• One such situation is regular-scale computing, when we perform computations on a usual computer, be it a single laptop of an affordable computer cluster.

• Regular-size computations work well in many cases, but, as we have mentioned earlier, there are many practical situations where regular-scale computing is not sufficient, where we need large-scale computing -what is usually called high-performance computing.

• For some situations, the opposite is true: we do not need even the full computational power of a regular laptop, it is sufficient to use limited computational power of an easierto-carry smaller device such as a cell phone or a smart watch. It is natural to call such situations small-scale computing.

Let us describe what are the typical resource limitations of these three types of situations -and what can we do to best takes these resource limitations into account.

Resource limitations of small-scale computing: energy. For small-scale computing -such as computing on cell phones -by definition, computational ability is not a problem. The main need for such small-scale devices comes from the fact that regular computers are not very portable: it is easier to carry a cell phone than a computer. This portability is a problem: to perform computations, we need energy. For a more or less stationary device like a regular computer, this is not a problem: we just plug in the computer, and get the energy non-stop. For portable devices, energy is a problem: we can only plug it in once in a while, and in a small volume, we can only a limited amount of energy. So, for such devices, the main resource limitation is energy.

This require a serious change in algorithms -since most traditional algorithms are designed to minimize computation time, without taking into account energy limitations. In particular, this means that energy-needing auxiliary procedures like garbage collection -that periodically frees the no-longer-used part of computer memory -have to be performed only based on need, and not, as in usual computers, periodically or whenever this leads to faster computations; see, e.g., [1] and references therein.

For example, a regular laptop usually updates when it is idle, when it is not doing any computations:

• this preserves its computational speed when it is used, and • no matter how many updates we make, it does not affect its computing abilities.

For portable devices, any such procedure uses some portion of stored energy, so it needs to be performed as rarely as possible.

The corresponding optimization of computations is a very challenging task, especially taking into account that some operating systems running on this devices are proprietary. Since the operating system uses a significant amount of energy, we thus face a partly-black-box problem even when the computational algorithms that we run on these devices are white-box ones; see, e.g., [2].

Resource limitations of regular-size computing: cost and cloud computing. For regularsize computing, computational ability is not a problem -otherwise, we would have needed a high-performance computer. As we have mentioned earlier, for such computers, energy is also not a problem: we just plug in. So is there any limited resource? Yes, a very mundane one: money.

The possibility to save money comes from the fact that the amount of needed computations changes with time:

• we have daily fluctuations, since most computations are performed at daytime,

• we have seasonal computations -e.g., a tax-advising company needs a lot of computations before the taxes deadline, stores need to perform more computations before Christmas, when many people are shopping for Christmas gifts, etc.

In all these cases, we have a choice:

• we can buy as many computers that we need at peak times, or

• we can buy fewer computers and at peak time, rent computation time on other computers.

Computations involving rented computation time are known as cloud computing; see, e.g., [3,4,5,6,7]. This name comes from the fact that we do not care where exactly these computations are performed -as long as they are performed. So, the possible locations of processors performing such computations do not form a clear shape, the resulting geographic shape is amorphous and often fast-changing, like a cloud in the sky. So how much computing power should we buy and how much should we rent? Detailed answers to this question can be found in [8,9,10,11]. In this paper, we show the optimal money-saving solution for the simplest case, when we have:

• full information about the costs and • full information about our computational needs. This means, in particular, that we know:

• the cost 𝑐 0 per computational step of in-house computations (including the corresponding part of the cost of buying the computers),

• the cost 𝑐 1 > 𝑐 0 per computational step of rented out-of-house computations, and

• probability distribution of the future computer needs 𝑥, which can be described, e.g., by the probability density function (pdf) 𝜌(𝑥).

Based on this information, we need to select the optimal amount 𝑥 0 of computational power that we should buy.

For each value 𝑥 > 𝑥 0 , we rent the missing amount 𝑥 − 𝑥 0 of the computational power. Thus, the overall expected cost of all these computations is equal to

𝑐 0 ⋅ 𝑥 0 + 𝑐 1 ⋅ ∫ ∞ 𝑥 0 𝑐 1 ⋅ (𝑥 − 𝑥 0 ) ⋅ 𝜌(𝑥) 𝑑𝑥.

To find the optimal value 𝑥 0 , we can differentiate this expression with respect to 𝑥 0 and equate the derivative to 0. As a result, we get the equation

𝐹 (𝑥 0 ) = 1 − 𝑐 0 𝑐 1 ,

where 𝐹 (𝑥) is the cumulative distribution function corresponding to the given pdf. This equation described the optimal amount of computing power that we should buy.

Comment. In [8,9,10,11], one can also find:

• algorithms for solving this problem in more realistic settings, when we know the costs and needs with some uncertainty,

• algorithms for deciding when a long-term contract is beneficial, and

• algorithms that help companies providing cloud computing to find the optimal locations of their computers around the world.

Resource limitations of large-scale computing: energy again. For large-scale computing, as we have mentioned earlier, one of the main problems is that some important practical problems still cannot be solved on such computers -computer engineers are working on this.

For the existing computers, what is the main limitation? Somewhat surprisingly, the main limitation is again energy, the same as for small-scale computing. Indeed, we can simply plug in a regular-size computer, but a usual company-wide or university-wide connection can only support a limited number of such plugged-in computers.

For a high-performance computer -which, crudely speaking, consists of hundreds and thousands of usual computers -we need to build a special power station, and even this may not be sufficient. How can we minimize the power needed for computations?

To answer this question, we need to go back to computer design:

• When we design a high-performance computer, we want to maximize the overall number of computations per second. From this viewpoint, a seemingly natural idea is to run all the processors at their maximum speed.

• It turns out that from the viewpoint of energy consumption, this is not the best arrangement.

To be more precise, if the power 𝑃 𝑖 used by each processor 𝑖 was exactly proportional to this processor's computation speed 𝑓 𝑖 -i.e., to the number of computational operations per second, this would not matter:

• if 𝑃 𝑖 = 𝑘 ⋅ 𝑓 𝑖 ,

• then, no matter how we distribute the task between processors, the overall power 𝑃 would be similarly proportional to the overall number 𝑓 of computational steps per second:

𝑃 = 𝑘 ⋅ 𝑓 .

In practice, however, the linear dependence of 𝑃 𝑖 on 𝑓 𝑖 is only a rough approximation, the actual dependence 𝑃 𝑖 = 𝐹 (𝑓 𝑖 ) is non-linear. As a result, instead of having all the processor work at their maximal speed, a better idea is:

• to have each processor work at the speed 𝑓 𝑖 at which the ratio

𝑃 𝑖 𝑓 𝑖 = 𝐹 (𝑓 𝑖 ) 𝑓 𝑖

is the smallest possible, and

• to use more processors if needed.

If at some point, the task requires a smaller overall computation speed 𝑓 , then, to minimize energy consumption:

• we make some processors idle, while

• others will work at their optimal speed (optimal in terms of energy consumption).

Quantum computing and its resource limitations. As we have mentioned several times, there is a practical need to increase computation speed beyond what we can do now. This need faces a fundamental physical limitation: that, according to relativity theory, the speed of all processes is limited by the speed of light. As a result:

• for a usual laptop of approximately 30 cm size, it takes 1 nanosecond for the light to travel from one side of the computer to another, and

• during this time, even the cheapest 4-GHz laptop performs 4 operations.

To drastically speed up processing, we need to drastically decrease the computer size -and this means to make memory cells drastically smaller. Already each cell consists of thousands of molecules. With further size decrease, we will get to the level of individual molecules -and at this level, we have to take into account physical laws of the micro-world -the law of quantum physics.

Computations that take quantum effects into account are know as quantum computing; see, e.g., [12]. One of the main features of quantum physics -see, e.g., [13,14] -is that, in addition to classical states 𝑠, 𝑠 ′ , . . . -which, in quantum physics, are denoted by |𝑠⟩, |𝑠 ′ ⟩, etc. -we can also have superpositions of these states, i.e., states of the type

𝑐 𝑠 |𝑠⟩ + 𝑐 𝑠 ′ |𝑠 ′ ⟩ + … ,

where 𝑐 𝑠 , 𝑐 𝑠 ′ , . . . , are complex values for which

|𝑐 𝑠 | 2 + |𝑐 𝑠 ′ | 2 + … = 1.

If we measure, in this superposition state, a classical variable 𝑣, then we get the value 𝑣(𝑠) with probability |𝑐 𝑠 | 2 , the value 𝑣(𝑠 ′ ) with probability |𝑐 𝑠 ′ | 2 , etc. The sum of all these probabilities should add to 1 -which explains the above condition on the values 𝑐 𝑠 , 𝑐 𝑠 ′ , . . . In particular, a quantum analogue of the usual bit -i.e., an element that can be in two possible states 0 and 1 -is a quantum bit (qubit, for short), i.e., the state of the type

𝑐 0 |0⟩ + 𝑐 1 |1⟩.

Similarly, a quantum analogue of a 2-bit state is a general 2-qubit state

𝑐 00 |00⟩ + 𝑐 01 |01⟩ + 𝑐 10 |10⟩ + 𝑐 11 |11⟩, etc.

There are many efficient quantum algorithms, e.g., Grover's algorithm [15,16,12].

• This algorithm enables to find an element with a given property in an unsorted 𝑛-element array in time proportional to √ 𝑛.

• On the other hand, for non-quantum algorithms, not to miss the desired element, we need to look at all 𝑛 elements, so we need at least 𝑛 computational steps, and 𝑛 ≫ √ 𝑛.

Grover's algorithm requires 𝑛-qubit states. At this moment, the number of available qubits is the main limitation to quantum computing, the main limited resource. What happens if we only have 𝑠 < 𝑛 qubits, i.e., if we can only implement 𝑠-qubit states? In this case, we can:

• divide the 𝑛-element array into 𝑛/𝑠 fragments of size 𝑠, and • apply Grover's algorithms to each fragment; see, e.g., [17].

For each fragment, Grover's algorithm requires √ 𝑠 steps, so the overall computation time will be 𝑛 𝑠

⋅ √ 𝑠 = 𝑛 √ 𝑠

. This amount:

• is still smaller than the time 𝑛 needed in the non-quantum case, and

• decreases as the resource bound 𝑠 increases (and reaches the Grover's value 𝑡 ∼ √ 𝑛 when we reach the case 𝑠 = 𝑛).

Black-Box Computing: Related Resource Limitations

General case. For black-box computing, when computations include calls to a "black-box" (proprietary or classified) program, a natural idea is to minimize the number of such calls:

• For proprietary programs, when we often need to pay for each class, this minimizes the overall costs.

• For classified programs, when each call is an extra vulnerability, this minimizes the risk.

• For programs that take a long time to compute, this also minimizes the overall computation time.

What can we compute in this manner? Usually, commercial software provides a turn-key solution to the corresponding problem, a solution that does not require additional programming. There exist commercial software for deep learning, for solving partial differential equations, etc. This software produces the result 𝑦 = 𝑓 (𝑥 1 , … , 𝑥 𝑛 ) of processing the inputs 𝑥 1 , … , 𝑥 𝑛 , but what is usually does not do is produce the accuracy of this result. To be more precise, most black-box packages compute the result implicitly assuming that the inputs are exactly known. In practice, the inputs come from measurements, and measurements are never absolutely accurate: the result x𝑖 of measuring the 𝑖-th input are, in general, different from the actual (unknown) value 𝑥 𝑖 of this input.

Thus, when we apply the black-box algorithm 𝑓 to the measurement results x𝑖 , the result ỹ of this computation will be, in general, different from the desired value 𝑦 = 𝑓 (𝑥 1 , … , 𝑥 𝑛 ):

ỹ = 𝑓 ( x1 , … , x𝑛 ) ≠ 𝑦 = 𝑓 (𝑥 1 , … , 𝑥 𝑛 ).

An important question is: how accurate is the computation result? In other words, how big can the difference Δ𝑦 def = ỹ = 𝑦 be? This question is very practically important. For example, if we are prospecting for oil and we learned that the given area contains ỹ = 200 million tons, then:

• if this amount is 200 ± 20, this is good news and we should start drilling, but • if it is 200 ± 300, then maybe there is no oil at all, and it is better to perform additional measurements before we spend a large amount of money on possibly useless drilling. • In some cases, we know the probability distribution of each measurement error Δ𝑥 𝑖 , and we know that these measurement errors are independent. In many cases, this distribution is normal (Gaussian) with mean 0 and known standard deviation 𝜎 𝑖 .

Let us formulate this

• In other cases, we do not have any information about the probabilities, we only know the upper bound Δ 𝑖 on the absolute value of the corresponding measurement error: |Δ𝑥 𝑖 | ≤ Δ 𝑖 ; see, e.g., [18].

How can we use this information to find bounds on Δ𝑥 𝑖 ?

Seemingly natural idea. If we know the values 𝑐 𝑖 , then, e.g., for the normal distribution, we can conclude that the linear combination • For 𝑐 𝑖 ≥ 0, the term 𝑐 𝑖 ⋅ Δ𝑥 𝑖 is increasing, so its largest value is attained when Δ𝑥 𝑖 attains its largest value Δ 𝑖 . The corresponding largest value of the 𝑖-th term is 𝑐 𝑖 ⋅ Δ 𝑖 .

• For 𝑐 𝑖 ≤ 0, the term 𝑐 𝑖 ⋅ Δ𝑥 𝑖 is decreasing, so its largest value is attained when Δ𝑥 𝑖 attains its smallest value −Δ 𝑖 . The corresponding largest value of the 𝑖-th term is −𝑐 𝑖 ⋅ Δ 𝑖 .

In both cases, the largest value of the 𝑖-th term is |𝑐 𝑖 | ⋅ Δ 𝑖 , so the largest value Δ of the sum Δ𝑦 is equal to

Δ = 𝑛 ∑ 𝑖=1 |𝑐 𝑖 | ⋅ Δ 𝑖 .

Straightforward approach and its limitations. If we knew the values 𝑐 𝑖 , then we could use the above formulas and find the desired characteristic 𝜎 or Δ of the approximation error Δ𝑦.

The algorithm 𝑓 (𝑥 1 , … , 𝑥 𝑛 ) is given as a black box, we do not know the corresponding expression and thus, we cannot simply differentiate this expression. However, what we can do is use numerical differentiation, according to which, for sufficiently small value ℎ 𝑖 , we have

𝑐 𝑖 = 𝜕𝑓 𝜕𝑥 𝑖 ≈ 𝑓 ( x1 , … , x𝑖 + ℎ 𝑖 , x𝑖+1 , … , x𝑛 ) − ỹ ℎ 𝑖 .

This way, we get the desired expressions for 𝜎 and Δ -but we need to make 𝑛 + 1 calls to 𝑓 : the first call to compute ỹ and then 𝑛 more calls to compute 𝑛 partial derivatives 𝑐 𝑖 . In many data processing algorithms, we use thousands of inputs, and each call to 𝑓 may require minutes or even hours, so this is not a very feasible idea. What can we do? Case of probability distributions: Monte-Carlo simulations. In situations when we know the probability distributions, the answer is well known: use Monte-Carlo simulations. Namely, for some integer 𝑁 , we:

• 𝑁 times simulate the variables 𝛿𝑥 (𝑘) 𝑖 , 𝑘 = 1, … , 𝑁 , whose distribution is the same as the distribution of the corresponding measurement errors Δ𝑥 𝑖 , and then

• compute the values

𝛿 (𝑘) = 𝑓 ( x1 , … , x𝑛 ) − 𝑓 ( x1 − 𝛿𝑥 (𝑘) 1 , … , x𝑛 − 𝛿𝑥 (𝑘) 𝑛 ).

By using the same Taylor expansion as before, we conclude that

𝛿 (𝑘) = 𝑛 ∑ 𝑖=1 𝑐 𝑖 ⋅ 𝛿𝑥 (𝑘) 𝑖 .

Since the values 𝛿𝑥 (𝑘) 𝑖 have the same distribution as the measurement errors Δ𝑥 𝑖 , we can conclude that the distribution of the values 𝛿 (𝑘) is exactly the same as the desired distribution of the approximation errors Δ𝑦.

In particular, for the normal distribution, we can estimate 𝜎 as

𝜎 ≈ √ 1 𝑁 ⋅ 𝑁 ∑ 𝑘=1 ( 𝛿 (𝑘) ) 2 .

Good news is that the accuracy of this estimation is proportional to 1 √ 𝑁 (see, e.g., [19]) and does not depend on the number of inputs 𝑛.

For example, if we want to find 𝜎 with accuracy 1 √ 𝑁 = 20%, we need to take 𝑁 = 25 iterations and thus, 25 calls for the program 𝑓 -which is much smaller than thousands of calls needed for straightforward estimation.

What if we only know the upper bounds?

At first glance, in this case the situation is hopeless: we do not know what distribution to use. However, interestingly, a Monte-Carlo approach is possible in this case as well; see, e.g., [20,21]. The idea is to use Cauchy distribution, i.e., the distribution with probability density function

𝜌(𝑥) = 1 𝜋 ⋅ Δ ⋅ 1 1 + ( 𝑥 Δ ) 2 .

The use of Cauchy distribution is based on the following property of this distribution:

• if independent random variables Δ𝑥 1 , . . . , Δ𝑥 𝑛 are Cauchy-distributed with parameters Δ 𝑖 ,

• then the linear combination Thus, to find Δ:

• we simulate 𝛿𝑥 (𝑘) 𝑖 to be Cauchy-distributed with parameter Δ 𝑖 ; • then the values

𝛿 (𝑘) = 𝑓 ( x1 , … , x𝑛 ) − 𝑓 ( x1 − 𝛿𝑥 (𝑘)

1 , … , x𝑛 − 𝛿𝑥 (𝑘) 𝑛 ) are Cauchy-distributed with the desired value of the parameter Δ. We can then find this value Δ by applying, e.g., the Maximum Likelihood method [19].

To solve this problem, researchers are working on designing robotic boats; see, e.g., [26] and references therein. The idea is to place this boat at the origin of the small river and let is flow along the river and measure the corresponding parameters. The problem is that measurements and transmission of the resulting data use some energy, and since boat is autonomous, it can only carry a limited amount of energy. So, the question is: how to get the desired map of the river by using the smallest possible number of measurements.

The main idea behind this minimization is straightforward:

• when we are passing through a part of the trajectory where the depth (or any other quantity of interest) is approximately constant, then there is no sense of measuring this quantity too frequently; but

• when we come to the part of the trajectory where the quantity starts changing, then we must measure it with the highest possible frequency.

Let us describe the corresponding problem in precise terms. We know reasonably well the trajectory followed by a spaceship or the trajectory followed by a floating robotic boat or by another similar device. Let 𝑥 be the parameter describing the device's location on this trajectory -e.g., the length of the path that it has already covered along this trajectory. At some values 𝑥, we perform some measurements. We may want to measure some quantity in its location -e.g., the depth at this location. We may want to measure some quantity that depends not only 𝑥, but also on some other parameters 𝑦 1 , …. For example, for the boat, we may want, by sending a signal in several directions, to measure the depths not only right beneath the boat, but also on other locations (𝑥, 𝑦 1 ) on the line orthogonal to the trajectory. In general, based on the measurements, we want to find the while dependence 𝑞(𝑥, 𝑦 1 , …) of the measured quantity on all these parameters, and we want to measure all these values with accuracy 𝜀. How can we reach this goal by using the smallest possible number of measurements?

A natural idea is to use a simple approximation

𝑞(𝑥, 𝑦 1 , …) = 𝐶 1 ⋅ 𝑞 1 (𝑥, 𝑦 1 , …) + … + 𝐶 𝑘 ⋅ 𝑞 𝑘 (𝑥, 𝑦 1 , …)

to describe the local behavior of the corresponding field. In this approximation:

• the functions 𝑞 1 (𝑥, 𝑦 1 , …), . . . , 𝑞 𝑘 (𝑥, 𝑦 1 , …) are fixed, and

• the parameters 𝐶 1 , … , 𝐶 𝑘 change from location to location.

In many cases, the field is smooth, so it can be locally approximated by linear or quadratic functions. For example, if we are interested in the depth 𝑑(𝑥) at different points along the trajectory, then we can use approximations 𝑑(𝑥) ≈ 𝐶 1 + 𝐶 2 ⋅ 𝑥 or take 𝑑(𝑥) ≈ 𝐶 1 + 𝐶 2 ⋅ 𝑥 + 𝐶 3 ⋅ 𝑥 2 .

• In the first case, 𝑘 = 2, 𝑞 1 (𝑥) = 1, and 𝑞 2 (𝑥) = 𝑥.

• In the second case, 𝑘 = 3, 𝑞 1 (𝑥) = 1, 𝑞 2 (𝑥) = 𝑥, and 𝑞 3 (𝑥) = 𝑥 2 .

How we proceed depends on:

• whether we know the probability distribution of the measurement error -e.g., Gaussian with standard deviation 𝜎 -or

• whether we only know the bound Δ on the measurement error.

In both cases, we first perform 𝑘 + 1 (or more) measurements corresponding to the values 𝑥 (𝑖) , 𝑦 (𝑖) 1 , . . . , where 𝑖 = 1, 2, … , 𝑘 + 1. Let 𝑞 (𝑖) be the corresponding measurement results.

In the probabilistic case, we form the resulting system of approximate equations:

𝑞 (1) ≈ 𝐶 1 ⋅ 𝑞 1 ( 𝑥 (1) , 𝑦 (1) 1 , … ) + … + 𝐶 𝑘 ⋅ 𝑞 𝑘 ( 𝑥 (1) , 𝑦 (1) 1 , … ) ;

… 𝑞 (𝑘+1) ≈ 𝐶 1 ⋅ 𝑞 1 ( 𝑥 (𝑘+1) , 𝑦 (𝑘+1) 1 , … ) + … + 𝐶 𝑘 ⋅ 𝑞 𝑘 ( 𝑥 (𝑘+1) , 𝑦 (𝑘+1) 1 , … ) .

We can use the usual Least Squares method (see, e.g., [19]) to find the estimates Ĉ1 , … , Ĉ𝑘 for the corresponding parameters 𝐶 1 , … , 𝐶 𝑘 . This method also provides us with the variance of these estimates, and with the covariances describing correlation between these estimates.

Based on these variances and covariances, for all future values 𝑥 and for all corresponding values 𝑦 1 , …, we can compute the standard deviation for the predicted value

𝑞(𝑥, 𝑦 1 , …) = Ĉ1 ⋅ 𝑞 1 (𝑥, 𝑦 1 , …) + … + Ĉ𝑘 ⋅ 𝑞 𝑘 (𝑥, 𝑦 1 , …).

As long as this value is smaller than or equal to the desired accuracy 𝜀, we do not perform any additional measurements. The next group of measurements is performed at the first future value 𝑥 at which the standard deviation becomes equal to 𝜀.

In the case of upper bounds, for each future value 𝑥 and for all possible values of the variables 𝑦 1 , …, we find the range [𝑞, 𝑞] of possible values 𝑞(𝑥, 𝑦 1 , …) by solving the following two optimization problems. To find the upper bound 𝑞, we maximize the expression 𝐶 1 ⋅ 𝑞 1 (𝑥, 𝑦 1 , …) + … + 𝐶 𝑘 ⋅ 𝑞 𝑘 (𝑥, 𝑦 1 , …) under the constraints 𝑞 (1) − Δ ≤ 𝐶 1 ⋅ 𝑞 1 ( 𝑥 (1) , 𝑦 (1) 1 , … ) + … + 𝐶 𝑘 ⋅ 𝑞 𝑘 ( 𝑥 (1) , 𝑦 (1) 1 , … ) ≤ 𝑞 (1) + Δ; … 𝑞 (𝑘+1) − Δ ≤ 𝐶 1 ⋅ 𝑞 1 ( 𝑥 (𝑘+1) , 𝑦 (𝑘+1) 1 , … ) + … + 𝐶 𝑘 ⋅ 𝑞 𝑘 ( 𝑥 (𝑘+1) , 𝑦 (𝑘+1) 1 , … ) ≤ 𝑞 (𝑘+1) + Δ. To find the lower bound 𝑞, we minimize the same expression under the same constraints. In both optimization problems, we optimize a function which is linear in terms of the unknowns 𝐶 1 , … , 𝐶 𝑘 under constraints which are also linear in terms of the unknowns. Such optimization problems are known as linear programming problems. There exist efficient algorithms for solving such problems; see, e.g., [27].

The next measurement is performed at a location for which the prediction accuracy is 𝜀, i.e., for which the interval [𝑞, 𝑞] has the form [ q − 𝜀, q + 𝜀] for some q, i.e., equivalently, for which

𝑞 − 𝑞 = 2𝜀.

Measuring combinations of quantities: general description and the case study of Covid-19 testing. When we measure, we are interested in some property. In some situations, only a few objects have the property of interest. In such situations, the possibility to measure combinations of quantities corresponding to individual objects can save on the number of measurements.

For example, if we are prospecting for oil, we save a lot of efforts if, instead of drilling at each possible location, we perform some measurements of the area as a whole and only drill in those areas which, according to these combined measurements, have oil.

Another example is Covid-19 testing. In the ideal world we should test everybody, but the number of available test kits is limited. If instead of testing all 𝑁 people from a given area, we apply the test to a mixture of samples from a group of several (𝑠 1 ) people, we can decrease the required number of tests kits -since:

• we can dismiss all the groups where no one was Covid-positive, and

• we need to continue testing only folks from the remaining groups; see, e.g., [28,29,30].

In [30], we considered the arrangement when after the group testing, we individually test everyone from still-suspicious groups. However, if there are many such folks, a reasonable idea is to combine these remaining folks into new groups of size 𝑠 2 < 𝑠 1 , and test each new group, etc., until, after 𝑛 such stages, on the following stage 𝑛 + 1, we test all remaining possibly-positive folks individually. Let us analyze what is the optimal approach to such multi-stage group testing.

Let 𝑝 denote the proportion of people in a given area that test Covid-positive. This number can be (and is) estimated based on the preliminary random testing. Thus, overall, we have 𝑁 ⋅ 𝑝 Covid-positive folks.

The value 𝑝 is small, so for sufficiently small group sizes 𝑠 𝑘 , the probability that we have two Covid-positive folks in each group is small. Since most of 𝑁 ⋅ 𝑝 Covid-positive folks belong to different groups, after 𝑘 stages, we will have 𝑁 ⋅ 𝑝 different groups of size 𝑠 𝑘 -and thus, we have 𝑁 ⋅ 𝑝 ⋅ 𝑠 𝑘 possibly-positive folks. . This is true for all 𝑘, so the value

𝑠 𝑘 𝑠 𝑘+1

is the same for all 𝑘. Let us denote this ratio by 𝑞. Thus, the values 𝑠 𝑘 form a geometric progression with common ratio 𝑞.

At the end, we check individually, so 𝑠 𝑛+1 = 1. Thus, we get 𝑠 𝑛 = 𝑞, 𝑠 𝑛−1 = 𝑞 2 , . . . , until we get 𝑠 1 = 𝑞 𝑛 . Thus, 𝑛 = log 𝑞 (𝑠 1 ) = ln(𝑠 1 ) ln(𝑞)

. For these values 𝑠 𝑘 , each term 𝑁 ⋅ 𝑝 ⋅ 𝑠 𝑘 𝑠 𝑘+1

in the expression for the overall number of tests is equal to 𝑁 ⋅ 𝑝 ⋅ 𝑞, and there are 𝑛 such terms -as many as stages in the testing procedure. Thus, the overall number of tests is equal to

𝑁 𝑠 1 + 𝑛 ⋅ 𝑁 ⋅ 𝑝 ⋅ 𝑞 = 𝑁 𝑠 1 + ln(𝑠 1 ) ln(𝑞) ⋅ 𝑁 ⋅ 𝑝 ⋅ 𝑞.

Minimizing this expression with respect to 𝑞 means minimizing the ratio 𝑞 ln(𝑞)

. Differentiating this expression with respect to 𝑞 and equating the derivative to 0, we conclude that ln(𝑞)−𝑞⋅ 1 𝑞 = 0, i.e., that ln(𝑞) = 1 and thus, 𝑞 = 𝑒. Multiplying both sides by 𝑠 2 1 , dividing both sides by 𝑁 ⋅ 𝑝 ⋅ 𝑒, and moving negative terms to the other side, we conclude that 𝑠 1 = 1 𝑝 ⋅ 𝑒 . Thus, we arrive at the following recommendation.

Optimal testing schedule. We first combine folks into groups of size 𝑠 1 = 1 𝑝 ⋅ 𝑒 , and test combined samples from each group. All the people from groups whose combined sample tested negative are Covid-negative and thus, do not need any additional testing.

The remaining folks are combined in groups of size 𝑠 2 = 𝑠 1 𝑒

. We test combined samples from each group. All the people from groups whose combined sample tested negative are Covid-negative and thus, do not need any additional testing.

The remaining folks are combined in groups of size 𝑠 3 = 𝑠 1 𝑒 2 , etc. At the end, we reach the stage at which 𝑠 𝑛+1 ≈ 1, i.e., where the remaining folks need to be individually tested.

How many tests do we need for the optimal testing schedule. Substituting the expression Is this method optimal? This is indeed the best we can do -at least asymptotically the best. Indeed, we need to find 𝑁 ⋅ 𝑝 elements out of 𝑁 . The number of such possible arrangement is

.

This expression differs from our number of tests by a factor 𝑒 ⋅ ln(2) ≈ 1.88; so, modulo a small multiplicative constant, our method is indeed asymptotically optimal.