{'fraction_non_alphanumeric': 0.06354161987332105, 'fraction_numerical': 0.014756749472171062, 'mean_word_length': 3.802394563693237, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimating Linear and Quadratic forms via Indirect Observations

In this paper, we further develop the approach, originating in [14 (arXiv:1311.6765),20 (arXiv:1604.02576)], to"computation-friendly"hypothesis testing and statistical estimation via Convex Programming. Specifically, we focus on estimating a linear or quadratic form of an unknown"signal,"known to belong to a given convex compact set, via noisy indirect observations of the signal. Most of the existing theoretical results on the subject deal with precisely stated statistical models and aim at designing statistical inferences and quantifying their performance in a closed analytic form. In contrast to this descriptive (and highly instructive) traditional framework, the approach we promote here can be qualified as operational -- the estimation routines and their risks are yielded by an efficient computation. All we know in advance is that under favorable circumstances to be specified below, the risk of the resulting estimate, whether high or low, is provably near-optimal under the circumstances. As a compensation for the lack of"explanatory power,"this approach is applicable to a much wider family of observation schemes than those where"closed form descriptive analysis"is possible. The paper is a follow-up to our paper [20 (arXiv:1604.02576)] dealing with hypothesis testing, in what follows, we apply the machinery developed in this reference to estimating linear and quadratic forms.


Introduction
This paper can be considered as a follow-up to the paper [27] dealing with hypothesis testing for simple families -families of distributions specified in terms of upper bounds on their momentgenerating functions. In what follows, we work with simple families of distributions, but our focus is on estimation of linear or quadratic forms of the unknown "signal" (partly) parameterizing the distribution in question. To give an impression of our approach and results, let us consider the sub-Gaussian case, where one is given a random observation ω drawn from a sub-Gaussian distribution P on R d : with sub-Gaussianity parameters µ ∈ R d , Θ = Θ T ∈ R d×d affinely parameterized by "signal" x ∈ R m . The goal is, given observation ω "stemming" from unknown signal x known to belong to a given convex compact set X ⊂ R m , to recover the value at x of a given linear form g(·) : R m → R. The estimate g we build is affine function of observation; the coefficients of the function, same as an upper bound on the -risk of the estimate on X 1 stem from an optimal solution to an explicit convex optimization problem and thus can be specified in a computationally efficient fashion. Moreover, under mild structural assumptions on the affine mapping x → (µ, Θ) the resulting estimate is provably nearoptimal in the minimax sense (see Section 4 for details). The latter statement is an extension of the fundamental result of D. Donoho [10] on near-optimality of affine recovery of a linear form of signal in Gaussian observation scheme.
This paper contributes to a long line of research on estimating linear (see, e.g., [36,24,16,35,29,26,7] and references therein) and quadratic ([21, 25, 1, 17, 14, 6, 15, 30, 18, 23, 32, 31, 28, 8] among others) functionals of parameters of probability distributions via observations drawn from these distributions. In the majority of cited papers, the objective is to provide "closed analytical form" lower risk bounds for problems at hand and upper risk bounds for the proposed estimates, in good cases matching the lower bounds. This paradigm can be referred to as "descriptive;" it relies upon analytical risk analysis and estimate design and possesses strong explanation power. It, however, imposes severe restrictions on the structure of the statistical model, restrictions making the estimation problem amenable to complete analytical treatment. There exists another, "operational," line of research, initiated by D. Donoho in [10]. The spirit of the operational approach is perfectly well illustrated by the main result of [10] stating that when recovering the linear form of unknown signal x known to belong to a given convex compact set X via indirect Gaussian observation ω = Ax + ξ, ξ ∼ N (0, I), the worstcase, over x ∈ X, risk of an affine in ω estimate yielded by optimal solution to an explicit convex optimization problem is within the factor 1.2 of the minimax optimal risk. Subsequent "operational" literature is of similar spirit: both the recommended estimate and its risk are given by an efficient computation (typically, stem from solutions to explicit convex optimization problems); in addition, in good situations we know in advance that the resulting risk, whether large or small, is nearly minimax optimal. The explanation power of operational results is almost nonexisting; as a compensation, the scope of operational results is usually much wider than the one of analytical results. For example, the just cited result of D. Donoho imposes no restrictions on A and X, except for convexity and compactness of X; in contrast, all known to us analytical results on the same problem subject (A, X) to severe structural restrictions. In terms of the outlined "descriptive -operational" dichotomy, our paper is operational. For instance, in the problem of estimating linear functional of signal x affinely parameterising the parameters µ, Θ of sub-Gaussian distribution we started with, we allow for quite general affine mapping x → (µ, Θ) and for general enough signal set X, the only restrictions on X being convexity and compactness.
Technically, the approach we use in this paper combines the machinery developed in [19,27] and the Cramer-type techniques for upper-bounding the risk of an affine estimate developed in [26]. 2 On the other hand, this approach can also be viewed as "computation-friendly" extension of theoretical results on "Cramer tests" supplied by [3,2,4,5] in conjunction with techniques of [13,14,11,12,10,8], which exploits the most attractive, in our opinion, feature of this line of research -potential applicability to a wide variety of observation schemes and (convex) signal sets X.
The rest of the paper is organized as follows. In Section 2 we, following [27], describe the families of distributions we are working with. We present the estimate construction and study its general properties in Section 3. Then in Section 4 we discuss applications to estimating linear forms of sub-Gaussian distributions. In Section 5 we apply the proposed construction to estimating quadratic forms of parameters of Gaussian and discrete distributions. To illustrate the performance of the proposed approach we describe results of some preliminary numerical experiments in which we compare the bounds on the risk of estimates supplied by our machinery with (numerically computed) lower bounds on the minimax risk. To streamline the presentation, all proofs are collected in the appendix.
Notation. In what follows, R n and S n stand for the spaces of real n-dimensional vectors and real symmetric n × n matrices, respectively; both spaces are equipped with the standard inner products, x T y, resp., Tr(XY ). Relation A B (A B) means that A, B are symmetric matrices of the same size such that A−B is positive semidefinite (resp., positive definite). We denote S n + = {S ∈ S n : S 0} and S n ++ = int S n + = {S ∈ S n : S 0}. We use "MATLAB notation:" [X 1 ; ...; X k ] means vertical concatenation of matrices X 1 , ..., X k of the same width, and [X 1 , ..., X k ] means horizontal concatenation of matrices X 1 , ..., X k of the same height. In particular, for reals x 1 , ..., x k , [x 1 ; ...; x k ] is a k-dimensional column vector with entries x 1 , ..., x k .
For probability distributions P 1 , ..., P K , P 1 × ... × P K is the product distribution on the direct product of the corresponding probability spaces; when P 1 = ... = P K , we denote P 1 × ... × P K by P K or [P ] K .
Given positive integer d, θ ∈ R d , Θ ∈ S d + , we denote by SG(θ, Θ) the family of all sub-Gaussian, with parameters (θ, Θ), probability distributions, that is, the family of all Borel probability distributions P on R d such that We use shorthand notation ω ∼ SG(θ, Θ) to express the fact that the probability distribution of random vector ω belongs to the family SG(θ, Θ).

Simple families of probability distributions
Let • F, 0 ∈ int F, be a closed convex set in Ω = R m symmetric w.r.t. the origin, • M be a closed convex set in some R n , • Φ(h; µ) : F × M → R be a continuous function convex in h ∈ F and concave in µ ∈ M.
Following [27], we refer to F, M, Φ(·, ·) satisfying the above restrictions as to regular data. Regular data F, M, Φ(·, ·) define the family of Borel probability distributions P on Ω such that We say that distributions satisfying (1) are simple. Given regular data F, M, Φ(·, ·), we refer to S[F, M, Φ] as to simple family of distributions associated with the data F, M, Φ. Standard examples of simple families are supplied by "good observation schemes," as defined in [26,19], and include the families of Gaussian, Poisson and discrete distributions. For other instructive examples and an algorithmic "calculus" of simple families, the reader is referred to [27]. We present here three examples of simple families which we use in the sequel.

Sub-Gaussian distributions
In this case, S[F, M, Φ] contains all sub-Gaussian distributions P on R m with sub-Gaussianity parameters from M: In particular, S[F, M, Φ] contains all Gaussian distributions N (θ, Θ) with (θ, Θ) ∈ M.

Quadratically lifted Gaussian observations
Let V be a nonempty convex compact subset of S d + . This set gives rise to the family P V of distributions of quadratic liftings [ζ; 1][ζ; 1] T of random vectors ζ ∼ N (θ, Θ) with θ ∈ R d and Θ ∈ V. Our goal now is to build regular data such that the associated simple family of distributions contains P V . To this end we select Θ * ∈ S d ++ and δ ≥ 0 such that for all Θ ∈ V one has Θ Θ * , and Θ 1/2 Θ Then (i) F, M + , Φ form a regular data, and for every (θ, Θ) ∈ R d × V it holds for all (h, H) ∈ F: (ii) Besides this, function Φ(h, H; Θ, Z) is coercive in the convex argument: whenever (Θ, Z) ∈ M + , For proof, see Appendix A.1.

Quadratically lifted discrete observations
Consider a random variable ζ ∈ R d taking values e i , i = 1, ..., d, where e i are standard basic orths in R d . 3 We identify the probability distribution P µ of such variable with a point µ = [µ 1 ; ...; Let now ζ K = (ζ 1 , ..., ζ K ) with ζ k drawn independently across k from P µ , and let We are about to point our regular data such that the associated simple family of distributions contains the distributions of the "quadratic lifts" ω[ζ K ] of random vectors ζ K .
and let Z d be a set of all positive semidefinite matrices from ∆ d . Denote Then for M = M (K) = K/2 , In other words, the simple family S[F, For proof, see Appendix A.2.

Situation and goal
Consider the situation as follows: given are Euclidean spaces E F , E M , E X along with • a nonempty set X contained in a convex compact set X ⊂ E X , • an affine mapping x → A(x) : E X → E M such that A(X ) ⊂ M, • a vector g ∈ E X and a constant c specifying the linear form G(x) = g, x + c : E X → R 4 , • a tolerance ∈ (0, 1).
Let P be the family of all Borel probability distributions on E F . Given a random observation where P ∈ P is associated with unknown signal x known to belong to X, "association" meaning that we want to recover the quantity G(x). Given ρ > 0, we call an estimate -a Borel function g(·) : E F → R -(ρ, )-accurate, if for all pairs x ∈ X, P ∈ P satisfying (11) it holds If ρ * is the infimum of those ρ for which estimate g is (ρ, )-accurate, then clearly g is (ρ * , )-accurate. We refer to ρ * as the -risk of the estimate g w.r.t. the data G(·), X, and (A, F, M, Φ): Risk ( g(·)|G, X, A, F, M, Φ) = min ρ : ∀(x, P ) ∈ X × P : When G, X, A, F, M, Φ are clear from the context, we shorten Risk ( g(·)|G, X, A, F, M, Φ) to Risk ( g(·)).
In the setting of this section, we are about to build, in a computationally efficient fashion, an affine estimate g(ω) = f * , ω + κ along with ρ * such that the estimate is (ρ * , )-accurate.

The construction
Let us set so that F + is a nonempty convex set in E F × R + , and let so that Ψ ± are convex real-valued functions on F + (recall that Φ is convex-concave and continuous on F × M, while A(X ) is a compact subset of M). These functions give rise to convex functions Ψ ± : E F → R given by and to convex optimization problem With our approach, a "presumably good" estimate of G(x) and its risk are given by an optimal (or nearly so) solution to the latter problem. The corresponding result is as follows: In the situation of Section 3.1, let Φ satisfy the relation Then and the functions Ψ ± (·) are convex real-valued. Furthermore, a feasible solutionf ,κ,ρ to the system of convex constraints in variables f , ρ, κ induces estimate of G(x), x ∈ X, with -risk at mostρ: Relation (17) (and thus -the risk bound (19)) clearly holds true whenf is a candidate solution to problem (13) andρ As a result, by properly selectingf we can make (an upper bound on) the -risk of estimate (18) arbitrarily close to Opt, and equal to Opt when optimization problem (13) is solvable.
For proof, see Appendix A.3.

Estimation from repeated observations
Assume that in the situation described in Section 3.1 we have access to K observations ω 1 , ..., ω K sampled, independently of each other, from a probability distribution P , and are allowed to build our estimate based on these K observations rather than on a single observation. We can immediately reduce this new situation to the previous one simply by redefining the data. Specifically, given It is immediately seen that the updated data satisfy all requirements imposed on the data in Section 3.1. Furthermore, for all f K = (f 1 , ..., f K ) ∈ F K , whenever a Borel probability distribution P on E F and x ∈ X satisfy (11), the distribution P K of K-element i.i.d. sample ω K = (ω 1 , ..., ω K ) drawn from P and x are linked by the relation Applying to our new data the construction from Section 3.2, we arrive at "repeated observations" version of Proposition 3.1. Note that the resulting convex constraints/objectives are symmetric w.r.t. permutations of the components f 1 , ..., f K of f K , implying that we lose nothing when restricting ourselves with collections f K with equal to each other components; it is convenient to denote the common value of these components f /K. With these observations, Proposition 3.1 becomes the statements as follows (we use the assumptions and the notation from the previous section): Proposition 3.2 In the situation described in Section 3.1, let Φ satisfy the relation (14), and let a positive integer K be given. Then functions Ψ ± : E F → R, are convex and real valued. Furthermore, letf ,κ,ρ be a feasible solution to the system of convex constraints in variables f , ρ, κ. Then, setting we get an estimate of G(x), x ∈ X, via independent K-repeated observations with -risk at mostρ, meaning that whenever a Borel probability distribution P is associated with x ∈ X in the sense of (11), one has Relation (21) clearly holds true whenf is a candidate solution to the convex optimization problem As a result, properly selectingf , we can make (an upper bound on) the -risk of estimate g(·) arbitrarily close to Opt, and equal to Opt when optimization problem (23) is solvable.
From now on, if otherwise is not explicitly stated, we deal with K-repeated observations; to get back to single-observation case, it suffices to set K = 1.

Application: estimating linear form of parameters of sub-Gaussian distributions 4.1 Situation
We are about to apply construction form Section 3 in the situation where our observation is sub-Gaussian with parameters affinely parameterized by signal x, and our goal is to recover a linear function of x. Specifically, consider the situation described in Section 3, with the data as follows: is the family of all sub-Gaussian distributions on R d ); • X = X ⊂ E X = R nx is a nonempty convex compact set, and Same as in Section 3, our goal is to recover the value of a given linear function G(y) = g T y + c at unknown signal x ∈ X via K-repeated observation ω K = (ω 1 , ..., ω K ) with ω i drawn, independently across i, from a distribution P which is associated with x, which now means "is sub-Gaussian with parameters (Ax + a, M (x))." We refer to Gaussian case as to the special case of the just described problem, where the distribution P associated with signal x is exactly N (Ax + a, M (x)).
In the case in question Φ(0; µ, M ) = 0, so that (14) takes place, and the left hand sides in the constraints (21) are Thus, system (21) reads We arrive at the following version of Proposition 3.2: Proposition 4.1 In the situation described above, given ∈ (0, 1), letf be a feasible solution to the convex optimization problem where Then the -risk of the affine estimate taken w.r.t. the data listed in the beginning of this section, is at mostρ.
It is immediately seen that optimization problem (24) is solvable, provided that and an optimal solution f * to the problem, taken along with yields the affine estimate with -risk, w.r.t. the data listed in the beginning of this section, at most Opt.

Consistency.
We can easily answer the natural question "when the proposed estimation scheme is consistent", meaning that for every ∈ (0, 1), it allows to achieve arbitrarily small -risk, provided that K is large enough. Specifically, if we denote G(x) = g T x + c, from Proposition 4.1 it is immediately seen that a sufficient condition for consistency is the existence off ∈ R d such thatf T Ax = g T x for all x ∈ X − X , or, equivalently, that g is orthogonal to the intersection of the kernel of A with the linear span of X − X . Indeed, under this assumption, for every fixed ∈ (0, 1) we clearly have lim K→∞ Φ(f ) = 0, implying that lim K→∞ Opt = 0, with Ψ and Opt given by (24). The condition in question is necessary for consistency as well, since when the condition is violated, we have Ax = Ax for properly selected x , x ∈ X with G(x ) = G(x ), making low risk recovery of G(x), x ∈ X , impossible already in the case of zero noise observations (i.e., those where the observation stemming from signal x ∈ X is identically equal to Ax + a) 5 .
Direct product case. Further simplifications are possible in the direct product case, where, in addition to what was assumed in the beginning of Section 4, It is immediately seen that in the direct product case problem (24) reads where Assuming v∈V Ker(M (v)) = {0}, the problem is solvable, and its optimal solution f * gives rise to the affine estimate Near-optimality. In addition to the assumption that we are in the direct product case, assume for the sake of simplicity, that M (v) 0 whenever v ∈ V . In this case (24) reads Now consider the problem of recovering g T u from observation ω i , 1 ≤ i ≤ K, independently of each other sampled from N (Au + a, M (v)), where unknown u is known to belong to U and v ∈ V is known. Let ρ (v) be the minimax -risk of the recovery: where inf is taken over all Borel functions g(·) : R Kd → R. Invoking [26, Proposition 4.1], it is immediately seen that whenever < 1 2 , one has where q N (s) is the s-quantile of the standard normal distribution. Since the family of all sub-Gaussian, with parameters (Au+a, we arrive at the following conclusion: Proposition 4.2 In the just described situation, the minimax optimal -risk Risk ( g(·)), random observations is within a moderate factor of the upper bound Opt on the -risk, taken w.r.t. the same data, of the affine estimate g * (·) yielded by an optimal solution to (26). Namely, with the "near-optimality factor"

Numerical illustration
In this section we consider the problem of estimating a linear form of signal x known to belong to a given convex compact subset X via indirect observations Ax affected by sub-Gaussian "relative noise." Specifically, our observation is Here A ∈ R d×n and Θ j ∈ S d + , j = 1, ..., n, are given matrices. In other words, we are in the situation where small signal results in low observation noise. The linear form to be recovered from observation The entities g, A, {Θ j } n j=1 and reals α ≥ 0 ("degree of smoothness"), σ > 0 ("noise intensity") are parameters of the estimation problem we intend to process. Parameters g, A, Θ j are generated as follows: • g ≥ 0 is selected at random and then normalized to have max x∈X g T x = 2; • we consider the case of n > d ("deficient observations"); the d nonzero singular values of A were set to where "condition number" θ ≥ 1 is a parameter; the orthonormal systems U and V of the first d left and, respectively, right singular vectors of A were drawn at random from rotationally invariant distributions; • positive semidefinite d × d matrices Θ j are orthogonal projectors on randomly selected subspaces in R d of dimension d/2 ; • in all experiments, we deal with single-observation case K = 1.
Note that X possesses ≥-largest pointx, whence M (x) M (x) whenever x ∈ X; as a result, sub-Gaussian distributions with matrix parameter M (x), x ∈ X, can be thought also to have matrix parameter M (x). One of the goals of the present experiment is to compare the risk of the affine estimate in the above model to its performance in the "envelope model" ω ∼ SG(Ax, M (x)), where the fact that small signals result in low-noise observations is ignored.
We present in Figure 1 the results of the experiment in which for a given set of parameters d, n, α, θ and σ we generate 100 random estimation problems -collections {g, A, Θ j , j ≤ d}. For each problem we compute (= 0.01)-risks of two affine in ω estimates of g T x as yielded by optimal solution to (24): the first -for the problem described above (the left boxplot in each group), and the second -for the aforementioned "direct product envelope" of the problem, where the mapping . Note the "noise amplification" effect (the risk is about 20 times the level σ of the observation noise) and significant variability of risk across the experiments. Seemingly, both these phenomena are due to deficient observation model (n > d) combined with "random interplay" between the directions of coordinate axes in R m (along these directions, X becomes more and more thin) and the orientation of the kernel of A.

Quadratic lifting and estimating quadratic forms
In this section we apply the approach in Section 3 to the situation where, given an i.i.d. sample ζ K = [ζ 1 ; ...; ζ K ], ζ i ∈ R d , with distribution P x of ζ i depending on an unknown "signal" x ∈ X, our goal is to estimate a quadratic functional q(x) = x T Qx + c T x of the signal. We consider two situations -the Gaussian case, where P x is a Gaussian distribution with parameters affinely depending on x, and discrete case where P x is a discrete distribution corresponding to the probabilistic vector Ax, A being a given stochastic matrix. Our estimation strategy is to apply the techniques developed in Section 3 to quadratic liftings ω of actual observations ζ (e.g., ω i = (ζ i , ζ i ζ T i ) in the Gaussian case), so that the resulting estimates are affine functions of ω's. We first focus on implementing this program in the Gaussian case.

Estimating quadratic forms, Gaussian case
In this section we focus on the problem as follows. Given are • a nonempty bounded set U ⊂ R m and a nonempty convex compact set V ⊂ R k , • an affine mapping u → A[u; 1] : R m → Ω = R d , where A is a given d × (m + 1) matrix, • a "functional of interest" where Q and q are known (m + 1) × (m + 1) symmetric matrix and k-dimensional vector, respectively.
We observe an i.i.d. sample ζ K = [ζ 1 ; ...; ζ K ], ζ i ∈ R d , with Gaussian distribution P u,v of ζ i depending on an unknown "signal" (u, v) known to belong to U × V : P u,v = N (A[u; 1], M (v)). Our goal is to estimate F (u, v) from observation ζ K . The -risk Risk ( g) of a candidate estimate g(·) -a Borel real-valued function on R Kd -is defined as the smallest ρ such that

Construction
Our course of actions is as follows.
• When adding to the above entities function Φ(·; ·)), as defined in (5), we conclude by Proposition 2.1 that M, F and Φ(·; ·) form a regular data such that for all (u, v) ∈ U × V and (h, H) ∈ F, where the inner product ·, · on E F is defined as (h, H), (g, G) = h T g + 1 2 Tr(HG), so that As a result of the above steps, we get at our disposal entities E X , E M , E F , F, M, Φ, X, X , A(·), G(·) and participating in the setup described in Section 3.1, and it is immediately seen that these entities meet all the requirements imposed by this setup. The bottom line is that the estimation problem stated in the beginning of this section reduces to the problem considered in Section 3.

The result
When applying to the resulting data Proposition 3.2 (which is legitimate, since Φ in (5) clearly satisfies (14)), we arrive at the result as follows: Proposition 5.1 In the just described situation, let us set (33) so that the functions Ψ ± (h, H) : R d × S d → R are convex. Furthermore, wheneverh,H,ρ,κ form a feasible solution to the system of convex constraints we get an estimate of the functional of interest with -risk not exceedingρ: In particular, setting for (h, we obtain an estimate (35) with -risk not exceedingρ.
For proof, see Section A.4.
Remark 5.1 In the situation described in the beginning of this section, let a set W ⊂ U ×V be given, and assume we are interested in recovering functional of interest (29) at points (u, v) ∈ W only. When reducing the "domain of interest" to W , we hopefully can reduce the -risk of recovery. Assuming that we can point out a convex compact set W ⊂ V × Z such that it can be straightforwardly verified that in this case the conclusion of Proposition 5.1 remains valid when the set V × Z in (33) is replaced with W, and the set U × V in (36) is replaced with W . This modification enlarges the feasible set of (34) and thus reduces the attainable risk bound.
Discussion. When estimating quadratic forms from K-repeated observations ζ K = [ζ 1 ; ...; ζ K ] with i.i.d. ζ i we applied "literally" the construction of Section 3.3, thus restricting ourselves with estimates affine in quadratic liftings ω i = (ζ i , ζ i ζ T i ) of ζ i 's. As an alternative to such "basic" approach, let us consider estimates which are affine in the "full" quadratic lifting ω = (ζ K , ζ K [ζ K ] T ) of ζ K , thus extending the family of candidate estimates (what is affine in ω 1 , ..., ω K , is affine in ω, but not vice versa, unless K = 1). Note that this alternative is covered by our approach -all we need, is to replace the original components d, M (·), V, A of the setup of this section with their extensions and set K to 1.
It is easily seen that such modification can only reduce the risk of the resulting estimates, the price being the increase in design dimension (and thus in computational complexity) of the optimization problems yielding the estimates. To illustrate the difference between two approaches, consider the situation (to be revisited in Section 5.2) where we are interested to recover the energy u T u of a signal u ∈ R m from observation where Θ is (unknown) diagonal matrix with diagonal entries from the range [0, σ 2 ], and a priori information about u is that u 2 ≤ R for some known R. Assume that (cf. Section 5.2.2) m ≥ 16 ln(2/ ), where ∈ (0, 1) is a given reliability tolerance and that R 2 ≥ mσ 2 . Under these assumptions one can easily verify that in the single-observation case the -risks of both the "plug-in" estimate ζ T ζ and of the estimate yielded by the proposed approach are, up to absolute constant factors, the same as the optimal -risk, namely, O(1)R, R = σ 2 m + σR ln(2/ ). Now let us look at the case K = 2 where we observe two independent copies, ζ 1 and ζ 2 , of observation (38). Here the -risks of the "naive" plug-in estimate 1 2 [ζ T i ζ 1 + ζ T 2 ζ 2 ], and of the estimate obtained by applying our "basic" approach with K = 2 are just by absolute constant factors better than in the single-observation case -both these risks still are O(1)R. In contrast to this, an "intelligent" plug-in 2-observation estimate ζ T 1 ζ 2 has risk O(1)σ(R + σ √ m) ln(2/ ) whenever R ≥ 0, which is much smaller than R when m ln(2/σ) and R ln(2/ ) σm. It is easily seen that with the outlined alternative implementation, our approach also results in estimate with "correct" -risk O(1)σ(R + σ √ m) ln(2/ ).

Consistency
We are about to present a simple sufficient condition for the estimator suggested by Proposition 5.1 to be consistent, in the sense of Section 4. Specifically, assume that A.1. V = {v} is a singleton such that M (v) 0, which allows to satisfy (3) with Θ * = M (v) and δ = 0, same as allows to assume w.l.o.g. that The consistency of our estimation procedure is given by the following simple statement: Proposition 5.2 In the just described situation and under assumptions A.1-2, given ∈ (0, 1), consider the estimate and Ψ ± = Ψ K ± are given by (33). Then the -risk of g K, (·) goes to 0 as K → ∞. For proof, see Section A.5.

The problem
Our first illustration is deliberately selected to be extremely simple: given direct noisy observation ζ = u + ξ of unknown signal u ∈ R m known to belong to a given set U , we want to recover the "energy" u T u of u; what we are interested in, is the quadratic in ζ estimate with as small -risk on U as possible; here ∈ (0, 1) is a given design parameter. Note that we are in the situation where the dimension d of the observation is equal to the dimension m of the signal underlying observation. The details of our setup are as follows: given. As a result, the "main ingredient" of constructions in Section 5.1 -the convex compact subset Z of Z + containing all matrices [u; 1][u; 1] T , u ∈ U , see (30), can be specified as • ξ ∼ N (0, Θ), with matrix Θ known to be diagonal with diagonal entries satisfying θσ 2 ≤ Θ ii ≤ σ 2 , 1 ≤ i ≤ d = m, with known θ ∈ [0, 1] and σ 2 > 0. In terms the setup of Section 5.1, we are in the case where V = {v ∈ R m : θσ 2 ≤ v i ≤ σ 2 , i ≤ m}, and M (v) = Diag{v 1 , ..., v m }; • the functional of interest is F (u, v) = u T u, i.e., is given by (29) with Q = I m and q = 0.

Processing the problem
It is easily seen that in the situation in question the construction in Section 5.1 boils down to the following: 1. We lose nothing when restricting ourselves with estimates of the form with properly selected scalars η and κ; 2. η and κ are supplied by the convex optimization problem (with just 3 variables α + , α − , η) where Specifically, the η-component of a feasible solution to (40) augmented by the quantity yields estimate (39) with -risk on U not exceeding Ψ(α + , α − , η); The "energy estimation" problem where ξ ∼ N (0, σ 2 I m ) with σ 2 known to belong to a given range is well studied in the literature. Available results investigate analytically the interplay between the dimension m of signal, the range of noise intensity σ 2 and the parameters R, r, and offer provably optimal, up to absolute constant factors, estimates. For example, consider the case with r = 0, and θ = 0, and assume for the sake of definiteness that R 2 ≥ σ 2 m (otherwise already the trivial -identically zero -estimate is near optimal) and that we are in "high dimensional regime," i.e., m ≥ 16 ln(2/ ). It is well known that in this case the optimal -risk, up to absolute constant factor, is σ 2 m + σR ln(2/ ) and is achieved, again, up to absolute constant factor, at the "plug-in" estimate x(ζ) = ζ T ζ. It is easily seen that under the circumstances similar risk bound holds true for the estimate (39) yielded by the optimal solution to (40). A nice property of the proposed approach is that (40) automatically takes care of the parameters and results in estimates with seemingly near-optimal performance, as is witnessed by the numerical results we present below.

Numerical results
In the experiments we are reporting on, we compute, for different sets of parameters m, r, R and θ (σ = 1 in all experiments) the 0.01-risk attainable by the proposed estimators in the Gaussian casethe optimal values of the problem (40), along with "suboptimality ratios" of such risks to the lower bounds on the best possible under circumstances 0.01-risks.
To compute these lower bounds we use the following construction. Consider the problem of estimating u 2 2 , u ∈ U = {u : r ≤ u 2 ≤ R} given observation ω = N (u, ϑI m ), with ϑ ∈ [θ, 1]. Same as in Section 4, the optimal -risk Risk opt for this problem is defined as the infimum of the -risk over all estimates. Now let us select somehow the r 1 , r 2 , r ≤ r 1 < r 2 ≤ R, and σ 1 , σ 2 , θ ≤ σ 1 , σ 2 ≤ 1, and let P 1 and P 2 be two distributions of observations as follows: P χ is the distribution of random vector ω = η + ξ, where η and ξ are independent, η is uniformly distributed over the sphere η 2 = r χ , and ξ ∼ N (0, σ 2 χ I m ), χ = 1, 2. It is immediately seen that if there is no test which can decide on the hypotheses H 1 : ω ∼ P 1 , and H 2 : ω ∼ P 2 via observation ω with total risk ≤ 2 (defined as the sum, over our two hypotheses, of probabilities to reject the hypothesis when it is true), the quantity is a lower bound on the optimal -risk Risk opt . In other words, denoting by p χ (·) the density of P χ , we have Now, the densities p χ are spherically symmetric, whence, denoting by q χ (·) the univariate density of the energy ω T ω of observation ω ∼ P χ , we have On a closest inspection, q χ is the convolution of two univariate densities representable by explicit computation-friendly formulas, implying that given r 1 , r 2 , σ 1 , σ 2 , we can check numerically whether the premise in (41) indeed takes place; whenever this is the case, the quantity is a lower bound on Risk opt 0.01 . In our experiments, we used a simple search strategy (not described here) aimed at crude maximizing this bound in r 1 , r 2 , σ 1 , σ 2 and used the resulting lower bounds on Risk opt 0.01 to compute the suboptimality ratios. 7 In Figures 2-4 we present some typical simulation results illustrating dependence of risks on problem dimension m (Figure 2), on ratio r/R (Figure 3), and on parameter θ (Figure 4). Different curves in each plot correspond to different values of the parameter R varying in {16, 32, 64, 128, 256, 512}, other parameters being fixed. We believe that quite moderate values of the optimality ratios presented in the figures (these results are typical for a much larger series of experiments we have conducted) attest a rather good performance of the proposed apparatus.

The problem
The estimation problem we address in this section is as follows. Our observations are where • P is a given d × m matrix, with m > d ("under-determined observations"), • u ∈ R m is a signal known to belong to a given compact set U , • ξ ∼ N (0, Θ) is the observation noise; Θ is positive semidefinite d × d matrix known to belong to a given convex compact set V ⊂ S d + . In our experiment, the data is specified as follows: 1. We assume that u ∈ R m is a discretization of a smooth function x(t) of continuous argument t ∈ [0; 1]:  3. The set V of allowed values of the "covariance" matrices Θ is the set of all diagonal d×d matrices with diagonal entries varying in [0, σ 2 ], with the "noise intensity" σ being a design parameter.

Processing the problem
Our estimating problem clearly is covered by the setups considered in Section 5.1. In terms of these setups, we specify Θ * as σ 2 I d , V as V, and M (v) as the identity mapping of S d onto itself; the mapping u → A[u; 1] becomes the mapping u → P u, while the set Z (which should be a convex compact subset of the set {Z ∈ S d+1 + : Z d+1,d+1 = 1} containing all matrices of the form [u; 1][u; 1] T , u ∈ U ) becomes the set As suggested by Proposition 5.1, linear in "lifted observation" ω = (ζ, ζζ T ) estimates of F (u) = 1 m u 2 2 stem from the optimal solution (h * , H * ) to the convex optimization problem with Ψ ± (·) given by (33) as applied with K = 1. The resulting estimate is and the -risk of the estimate is (upper-bounded by) Opt. Problem (43) is a well-structured convex-concave saddle point problem and as such is beyond the "immediate scope" of the standard Convex Programming software toolboxes primarily aimed at solving well-structured convex minimization problems. However, applying conic duality, one can easily eliminate in (33) the inner maxima over v, Z ro arrive at the reformulation which can be solved numerically by CVX [20], and this is how (43) was processed in our experiments.

Numerical results
To quantify the performance of the proposed approach, we present, along with the upper risk bounds, simple lower bounds on the best -risk achievable under the circumstances. The origin of these lower bounds is as follows. Let w ∈ U with t(w) = P w 2 , and let ρ = 2σq N (1 − ) where q N (·) is the standard normal quantile: Then for θ(w) = max[1 − ρ/t(w), 0], we have w := θ(w)w ∈ U , and P w − P w 2 ≤ ρ. The latter, due to the origin of ρ, implies that there is no test which decides on the hypotheses u = w and u = w via observation P u + ξ, ξ ∼ N (0, σ 2 I d ), with risk < . As an immediate consequence, the quantity is a lower bound on the -risk, on U , of a whatever estimate of u 2 2 . We can now try to maximize the resulting lower risk bound over U , thus arriving at the lower bound On a closest inspection, the latter problem is not a convex one, which does not prevent us from building its suboptimal solution.
Note that in our experiments even with fixed design parameters d, m, θ, σ, we still deal with families of estimation problems differing from each other by their "sensing matrices" P ; orientation of the system of right singular vectors of P with respect to the axes of U is random, so that these matrices varies essentially from simulation to simulation, which affects significantly the attainable estimation risks. We display in Figure 5 typical results of our experiments. We see that the (theoretical upper bounds on the) -risks of our estimates, while varying significantly with the parameters of the experiment, all the time stay within a moderate factor from the lower risk bounds.

Estimation of quadratic functionals of a discrete distribution
In this section we consider the situation as follows: we are given an d × m "sensing matrix" A which is stochastic -with columns belonging to the probabilistic simplex ∆ d = {v ∈ R d : v ≥ 0, i v i = 1}, and a nonempty closed subset U of ∆ m , along with a K-repeated observation ζ K = (ζ 1 , ..., ζ K ) with ζ i , 1 ≤ i ≤ K, drawn independently across i from the discrete distribution µ = Au * , where u * is an unknown probabilistic vector ("signal") known to belong to U . We always assume that K ≥ 2. We treat a discrete distribution on d-point set as a distribution P µ on the d vertices e 1 , ..., e d of ∆ d , so that possible values of ζ i are basic orths e 1 , ..., e d in R d with Prob ζ∼µ (ζ = e j ) = µ j . Our goal is to recover from observation ζ K the value at u * of a given quadratic form

Construction
Observe that for u ∈ ∆ m , we have u = [uu T ]1 m , where 1 m is the all-ones vector in R m . This observation allows to rewrite F (u) as a homogeneous quadratic form: Our goal is to construct an estimate g(ζ K ) of F (u), specifically, estimate of the form where ω[ζ K ] is the "quadratic lifting" of observation ζ K (cf. (6)): and h ∈ S m and κ ∈ R are the parameters of the estimate. To this end • we set x(u) = uu T , with X = {uu T : u ∈ U }, and specify a convex compact subset X of the intersection of the "symmetric matrix simplex" ∆ m ⊂ S m (see (7)) and the cone S m + of positive semidefinite matrices such that X ⊂ X ⊂ E X := S m . We put F = E F := S d , and M = ∆ d , thus AX A T ⊂ M ⊂ E M := S d .
• By Proposition 2.2, F, M and Φ(·; ·), as defined in (8), form a regular data such that setting M = K/2 , for all u ∈ U and h ∈ S d it holds where h, w = Tr(hw) is the Frobenius inner product on S d .
Observe that for x ∈ E X , x → A(x) = AxA T is an affine mapping from X into M, and setting we get a linear functional on E X such that we ensure that The relation Φ(0, z) = 0 ∀z ∈ M being obvious, Proposition 2.2 combines with Proposition 3.1 to yield the following result. The functions Ψ ± are real valued and convex on S m , and every candidate solutionh to the convex optimization problem induces the estimate of the functional of interest (45) via observation ζ K with -risk on U not exceedingρ = Ψ(h):

Numerical illustration
To illustrate the above construction, consider the following problem: we observe independent across k ≤ K realizations ζ k of discrete random variable ζ taking values 1, ..., d. The distribution p ∈ ∆ d of ζ is linearly parameterized by "signal" u which itself is a probability distribution on "discrete square" Ω = Ξ × Ξ, Ξ = {1, ..., m}: Here A i,rs ≥ 0 are known coefficients such that i A i,rs = 1 for all (r, s) ∈ Ω. Now, given two sets I ⊂ Ξ and J ⊂ Ξ, consider the events I = I × Ξ ⊂ Ω and J = Ξ × J ⊂ Ω. Our objective is to quantify the deviation of these events, the probability distribution on Ω being u, from independence, specifically, to estimate, via observations ζ 1 , ..., ζ K , the quantity which is a quadratic function of u. In the experiments we report below, this estimation was carried out via a straightforward implementation of the construction presented earlier in this section. Our setup was as follows: 1. We use d = m 2 . d×d column-stochastic "sensing matrix" A 8 corresponding to the "mixed-noise observations" [33,34] is generated according to A = θI d + (1 − θ)D, with column-stochastic d × d matrix D, θ ∈ [0, 1] being our control parameter. D was selected at random, by normalizing columns of a d×d matrix with independent entries drawn from the uniform distribution on [0, 1]; 2. We set X = {x ∈ S d : x rs,r s ≥ 0 ∀r, s, r , s ≤ m, x 0, We present in Figure 6 the results of experiments for θ taking values in {0.00, 0.25, 0.50, 0, 75, 1.00}. Other things being equal, the smaller θ, the larger is the condition number cond(A) of the sensing matrix, and thus the larger is the (upper bound on the) risk of our estimate -the optimal value of (47). Note that the variation of F ij over X is exactly 1/2, so the maximal risk is ≤ 1/4. It is worthy to note that simple (if compared, e.g., to much more involved results of [22]) bounds in Proposition 2.2 for Laplace functional of order-2 U -statistics distribution result in fairly good approximations of the risk of our estimate (cf. the boxplots of empirical distributions of the estimation error in the right plot of Figure 6).
Besides this, setting F (X) = − ln Det(X) : int S d + → R and equipping S d with the Frobenius inner product, we have ∇F (X) = −X −1 , so that with R 0 = Θ 1/2 * HΘ 1/2 * , R 1 = Θ 1/2 HΘ 1/2 , and ∆ = R 1 −R 0 , we have for properly selected λ ∈ (0, 1) and R λ = λR 0 + (1 − λ)R 1 : We conclude that Denoting by µ i the eigenvalues of R λ and noting that R λ ≤ max[ R 0 , R 1 ] = d(H) we get |µ i | ≤ d(H). Therefore, the eigenvalues (52), we conclude that and when substituting into (53) we get Furthermore, because by (3) Consequently,   To complete the proof of (i), all we need is to verify the claim that F, M + , Φ is regular data, which boils down to checking that Φ : F × M + → R is continuous and convex-concave. Let us verify convexity-concavity and continuity. Recalling that Υ(H; Θ) : F × V → R indeed is convex-concave and continuous, the verification in question reduces to checking that Γ(h, H; Z) is convex-concave and continuous on (R d × H γ ) × Z + . Continuity and concavity in Z being evident, all we need to prove is that whenever Z ∈ Z + , the function Γ Z (h, H) := Γ(h, H; Z) is convex in (h, H) ∈ F = R d × H γ . By the Schur Complement Lemma, we have where N is the number of permutations σ ∈ S K such that a particular pair (i, j), 1 ≤ i = j ≤ K is met among the pairs (σ 2k−1 , σ 2k ), 1 ≤ k ≤ M . Comparing the total number of ω ij -terms in the left and the right hand sides of the latter equality, we get Card(S K )M = N K(K − 1), which combines with the equality itself to imply that Let σ id be the identity permutation of 1, ..., K. Due to (55) we have The distribution of the random variable ω 12 [ζ 2 ] = 1 2 [ζ 1 ζ T 2 + ζ 2 ζ T 1 ], ζ 2 ∼ P 2 µ , clearly is P Z[µ] , so that ln E ζ 2 ∼P The latter relation combines with (56) to imply (9).