Sparse Quadratic Logistic Regression in Sub-quadratic Time

TL;DR

This paper introduces a fast, sub-quadratic algorithm for support recovery in sparse quadratic logistic regression, leveraging novel correlation tests and hashing techniques to efficiently identify relevant features.

Contribution

The paper presents a new algorithm that significantly reduces computation time for support recovery in sparse quadratic logistic regression by combining correlation tests and hashing methods.

Findings

Algorithm achieves sub-quadratic time complexity.

Effective support recovery demonstrated through experiments.

Novel correlation test for non-binary finite support cases.

Abstract

We consider support recovery in the quadratic logistic regression setting - where the target depends on both p linear terms $x_i$ and up to $p^2$ quadratic terms $x_i x_j$. Quadratic terms enable prediction/modeling of higher-order effects between features and the target, but when incorporated naively may involve solving a very large regression problem. We consider the sparse case, where at most $s$ terms (linear or quadratic) are non-zero, and provide a new faster algorithm. It involves (a) identifying the weak support (i.e. all relevant variables) and (b) standard logistic regression optimization only on these chosen variables. The first step relies on a novel insight about correlation tests in the presence of non-linearity, and takes $O(pn)$ time for $n$ samples - giving potentially huge computational gains over the naive approach. Motivated by insights from the boolean case, we propose a non-linear correlation test for non-binary finite support case that involves hashing a variable and then correlating with the output variable. We also provide experimental results to demonstrate the effectiveness of our methods.

Figures

Tables

Table 1. Table 1: Comparison of standard classifiers on the Dorothea dataset with and without addition of interaction terms using our correlation test: We take two standard classifiers -Linear SVM and Logistic Regression both with ℓ 1 subscript ℓ 1 \ell_{1} regularization. We train them on the binary feature matrix as is. Subsequently, we use the correlation test from Algorithm 1 to identify weak support (controlled by a sparsity parameter k 𝑘 k ). We create new features by multiplying binary features in the weak support pairwise and train the classifiers again to compare.

Classifier	AUC Score	Sparsity after CV	Running Time
Logistic + ${\ell }_1$	0.890357	N/A	6.197472s
Logistic + ${\ell }_1$ + quad terms	0.909252	8	12.022856s
LSVM + ${\ell }_1$	0.822040	N/A	3.986561s
LSVM+ ${\ell }_1$ +quad terms	0.881608	10	32.880251s

Table 2. Table 2: Average fraction of the weak support recovered vs samples for the non-linear correlation test. Nonlinear correlation of each variable is the average of the correlation obtained using 10 random hashes. Weak support recover rate is the fraction of the weak support that is in the top 20 variables with highest nonlinear correlation.

Samples	$500$	$1000$	$5000$	$10000$
WS
Rec. Rate	$0.346$	$0.488$	$0.732$	$0.822$

Equations

$$\mathrm{Pr}(Y=1|\mathbf{X}=\mathbf{x})=\sigma (\gamma f(\mathbf{x}))=\frac{\exp (\gamma f(\mathbf{x}))}{1+\exp (\gamma f(\mathbf{x}))},$$

$$\text{where }f(x)=\sum _{(i,j)\in Q}{\beta }_{i,j}{x}_i{x}_j+\sum _{j\in L}{\alpha }_j{x}_j+c$$

$$\mathbb{E}\left[Y\frac{(Z-E[Z])}{\mathrm{stddev}(Z)}\right]>\theta .$$

$$\sigma (x)$$

$$L(x)$$

$$\mathbb{E}[Y(X_k-{\mu }_k)]$$

$$=\mathbb{E}[(X_k-{\mu }_k)\mathbb{E}[Y|\mathbf{X}]]$$

$$=\mathbb{E}[\sigma (\gamma f(\mathbf{X}))(X_k-{\mu }_k)]$$

$$\mathbb{E}[Y(X_k-{\mu }_k)]=\sum _{i=1}^{m+1}\sigma (\gamma v_i)[1_{i\le m}+1_{i=m+1,E_0=1}]$$

$$\dots [p_k(1-{\mu }_k)\mathrm{Pr}(\mathbf{x}:f(\mathbf{x})=v_i|x_k=+1)$$

$$\dots -(1+{\mu }_k)(1-p_k)\mathrm{Pr}(\mathbf{x}:f(\mathbf{x})=v_i|x_k=-1)]$$

$$=2p_k(1-p_k)\left[\sigma (0)g(0)1_{E_0=1}+\sum _{i=1}^m\sigma (\gamma v_i)g(v_i)\right]$$

$$\mathrm{Pr}\left(\sum _{k\in [n]}{\alpha }_i{b}_i\in (-\nu /2,\nu /2)\right)\le \frac{\left(\genfrac{}{}{0ex}{}{n}{\lfloor n/2\rfloor }\right)}{2^n}\approx \frac{1}{\sqrt{n}}$$

$$\sigma (0)g(0)1_{E_0=1}+\sum _{i=1}^m\sigma (\gamma v_i)g(v_i)=d_0[\sum _{i=1}^{m}g(v_i)+$$

$$g(0)1_{E_0=1}]+\sum _{j=0}^{\infty }{d}_{2j+1}{\gamma }^{2j+1}\left[\sum _{i=1}^m{v}_i^{2j+1}g(v_i)\right]$$

$$\displaystyle\left[\begin{array}[]{cccc}v_{1}&v_{2}&v_{3}&\ldots v_{m}\\ v_{1}^{3}&v_{2}^{3}&v_{3}^{3}&\ldots v_{m}^{3}\\ \cdots&\cdots&\cdots&\cdots\\ v_{1}^{2i+1}&v_{2}^{2i+1}&v_{3}^{2i+1}&\ldots v_{m}^{2i+1}\\ \cdots&\cdots&\cdots&\cdots\\ \end{array}\right]\left[\begin{array}[]{c}g(v_{1})\\ g(v_{2})\\ \cdots\\ g(v_{m})\\ \end{array}\right]=0$$

$$\exists v_i:g(v_i)\ne 0\mathrm{and}\mathrm{for}\mathrm{the}\mathrm{same}{v}_i\mathrm{if}\exists j:v_j=-v_i,$$

$$\mathrm{then}g(v_j)\ne g(v_i).$$

$$g(v_i)$$

$$x_j=y_j\forall j\in L)\ast [\mathrm{Pr}\left(x_{\ell }=z_{\ell }^{+1}\forall \ell \in C-\{1\}\right)$$

$$-\mathrm{Pr}\left(x_{\ell }=z_{\ell }^{-1}\forall \ell \in C-\{1\}\right)]$$

$$=\mathrm{Pr}(x_i{x}_j=w_{i,j}\forall (i,j)\in Q:i,j\notin C,$$

$$x_j=y_j\forall j\in L)\ast \left(\prod _{j\in C-\{1\}}{p}_j\right)\ast$$

$$\left[\prod _{j\in C-\{1\}}{\left(\frac{1-p_j}{p_j}\right)}^{1-q_j^{+}}-\prod _{j\in C-\{1\}}{\left(\frac{1-p_j}{p_j}\right)}^{q_j^{+}}\right]$$

$$e(v_i)$$

$$x_j=y_j\forall j\in L)$$

$$|g(v_i)|$$

$$|\log (1/p_i-1)|>{\delta }^{\prime }.$$

$$|g(v_i)|$$

$$\left|\prod _{j\in C-\{1\}}{\left(\frac{1-p_j}{p_j}\right)}^{1-2q_j^{+}}-1\right|$$

$$\ge {\delta }^r{\left(\frac{1}{1-\delta }-1\right)}^{|C|-1}\left|e^{\sum _j\log \left(\frac{1-p_j}{p_j}\right)\left(1-2q_j^{+}\right)}-1\right|$$

$$\ge {\delta }^{2r}\left|e^{\sum _{j\in C-\{1\}}\log \left(\frac{1-p_j}{p_j}\right)\left(1-2q_j^{+}\right)}-1\right|$$

$$\mathrm{Pr}(|g(V)|>{\delta }^{2r}\min (2^{{\delta }^{\prime }/2}-1,1-2^{-{\delta }^{\prime }/2})$$

$$\text{hfill}\ge 1/2.$$

$$\mathbb{E}[Y(X_k-{\mu }_k)]$$

$$\sum _{i=1}^m\sigma (\gamma v_i)g(v_i)]$$

$$=p_1(1-p_1)\left[1_{E_0=1}g(0)+\sum _{i=1}^{m}g(v_i)\right]$$

$$+p_1(1-p_1)\left[\sum _{i=1}^{m}L(\gamma v_i)g(v_i)\right]$$

$$\mathbb{E}[Y(X_k-{\mu }_k)]$$

$$=p_1(1-p_1)\left[\sum _{i=1}^{m}L(\gamma |v_i|)g^{\prime }(v_i)\right]$$

$$\left[\begin{array}[]{cccc}|v_{1}|&|v_{2}|&\cdots&|v_{m}|\\ 1&|v_{2}|&\cdots&|v_{m}|\\ 1&1&\cdots&|v_{m}|\\ \cdots&\cdots&\cdots&\cdots\\ 1&1&\cdots&1\end{array}\right]\left[\begin{array}[]{c}g^{\prime}(v_{1})\\ \cdots\\ \cdots\\ g^{\prime}(v_{m})\end{array}\right]=\left[\begin{array}[]{c}b_{0}\\ \cdot\\ \cdot\\ b_{m}\end{array}\right]$$

$$X_1:$$

Full text

1 Introduction

In this paper we consider the following simple quadratic logistic regression setting: we have $p$ predictor variables $X_{i}\in\mathcal{X}\subset\mathbb{R}$ each with finite support on the real line($\lvert\mathcal{X}\rvert$ is finite), with $i=1,\ldots,p$, and one binary target variable $Y\in\{0,1\}$ which depends on $X$ as follows:

[TABLE]

Note that, here we have indexed the quadratic terms by $Q\subseteq[p]\times[p]$ and the linear terms by $L\subseteq[p]$. Here, $\sigma(\cdot)$ is a monotone non-linear function whose output range is in $[0,1]$ and $\gamma$ is a rescaling constant. As a specific example, we can consider the non-linear function $\sigma$ to be a logistic function $\sigma(\cdot)=\frac{\exp(\cdot)}{1+\exp(\cdot)}$, but our results hold more generally. Here, $\gamma\in\mathbb{R}^{+}$ is a scaling parameter.

Let the coefficients be bounded, i.e., $0<a\leq\lvert\beta_{i,j}\rvert,~{}\lvert\alpha_{j}\rvert,~{}\lvert c\rvert\leq b$. We assume that the function $f$ is $s$-sparse, i.e. $\lvert Q\rvert+\lvert L\rvert\leq s$. The number of variables that the function depends on is $s\leq r\leq 2s$. This problem is called a sparse logistic regression with real features with finite support. For example, $y$ can indicate if a user will click on a displayed advertisement and $x_{i}$ can be known user features. The function $f(\mathbf{x})$ captures how desirable an ad is to a user (each ad will be described by its own polynomial coefficients $\alpha,\beta$) and the non-linear function translates that into a probability of clicking which is then observed through a Bernoulli random variable $y$. The problem is to learn the polynomial coefficients $\alpha,\beta$ from samples $(\mathbf{x},y)$ generated from some known distribution.

The well-known challenge with learning higher-order polynomials is the sample and the computational complexity. The model has $\Theta(p^{2})$ coefficients that must be learned, and the number of features $p$ can be in the order of thousands or even millions for many applications. The key simplifying assumption is that the number of actually influential features (and influential pairs of features) can be small. We therefore assume that the total number of non-zero linear and bi-linear terms in $f(\mathbf{x})$ is at most $s$, i.e. $||\alpha||_{0}+||\beta||_{0}\leq s$.

Several methods have been used for sparse logistic regression, including maximum likelihood estimation with regularization and greedy methods [15, 10, 11]. All previously known approaches for this problem have a running time that is $\Theta(p^{2}n)$, i.e. quadratic in the ambient dimension $p$.

1.1 Our Contributions

Algorithm: We propose two algorithms for sparse quadratic support recovery for logistic regression, one for binary input variables and the other for bounded non-binary real valued variables. The algorithm for the binary case is (a) a simple correlation test per variable that identifies relevant variables (termed the weak support) and (b) a standard ML optimization on the set of all linear and quadratic terms formed by the weak support. The first step takes $O(pn)$ time for $p$ variables and $n$ samples; the second step is standard logistic regression with $O(w^{2})$ features, where $w$ is the size of the weak support. The main insight here is that except for a measure-zero set of coefficients, standard correlation tests recover the weak support even for quadratic terms – but only when there are non-linearities (like the sigmoid function in logistic regression). For the case of non-binary real valued variables with finite support, we propose a generalized non-linear correlation test for weak support recovery. First, a candidate variable is hashed using a random hash function and a correlation test performed to recover the weak support. The runtime complexity of the first step in this case is also $O(pn)$.

Analysis: In the binary case, we show that a simple linear correlation test is asymptotically effective when the non-linear function $\sigma(\cdot)$ is strongly non-linear, i.e. it has $2^{s}$ non-zero terms in its Taylor expansion, and the polynomial coefficients are not in a set of measure 0 which we precisely characterize. We show finite sample complexity results for a specific non-linearity (piecewise-linear $\sigma(\cdot)$). In short, the fundamental property that is key to our results is the non-linearity $\sigma(\cdot)$ and an extra degree of freedom $\gamma$. Insights about non-linearity plays a strong role behind our non-linear correlation test for real variables with finite support. This test introduces non-linearity by transforming the features using a hash function before calculating the correlation with the target variable.

Experiments: We show that linear correlation test works on synthetic randomly chosen functions with tens of relevant variables and $1000$s of irrelevant variables. With tens of thousands of samples, it recovers $>90\%$ of the weak support. We also show how correlation test and addition of few quadratic terms based on the weak support to the set of features helps improve standard classifiers on a real data set. We also similar positive experimental results for the non-linear correlation test in the non-binary case.

1.2 Related Work

Logistic regression is a very well studied problem in statistics. Logistic regression involving a linear function of the variables is carried out by computing a maximum-likelihood estimator which is known to be convex. If the number of variables is small, i.e., if we have a sparse logistic regression problem, then $\ell_{1}$ regularization (sometimes other regularizers) is used to enforce sparsity. Obtaining faster convex optimization based solvers for this particular convex optimization problem is the focus for many prior studies (see [11][10] and [6]). In the online setting, [1] consider the problem of adaptively including higher order monomials when using online convex optimization tools for regression problems on polynomials.

Another important classic problem is learning Ising models with pairwise interactions. Learning pairwise Ising Models on $p$ nodes can be thought of as $p$ parallel logistic regression problems involving a linear functional. Ravikumar et al. [16] propose convex optimization methods to solve these problems while Bresler et al. [2] showed that simple correlation tests and a greedy procedure is sufficient. The second work is closer to ours in spirit. However, we work with quadratic polynomials of independent variables with finite support.

Another line of work in statistics literature [13, 4] investigates the problem of identifying relevant linear and quadratic terms in a logistic regression problem in high dimensional settings with real covariates. [4] propose a screening procedure to identify relevant variables. However, they consider a Gaussian mixture model and their method needs estimation of precision matrices that has quadratic runtime complexity. The recent work of [13] proposes a greedy selection procedure for interaction terms and linear terms. At every step, their method involves evaluating the change in the logistic likelihood by adding every possible variable outside the current set and all the relevant quadratic terms. This has linear complexity per step. The evaluation of logistic likelihood has to be done on a growing set of variables which is also expensive. Contrary to these prior works, in this paper, we focus on real variables with finite support and achieve a total running time that is linear in the number of variables.

There has been a lot of work analyzing the time and sample complexity requirements for learning sparse polynomials over the boolean domain when the evaluations of the polynomial over the samples are directly observed. The closest to our work are [12, 9]. In [9], the authors study the general problem of learning a sparse polynomial function from random boolean samples. They define the unique sign property (USP) property which is also key for our machinery. However, it has several key differences since it does not achieve sub-quadratic learning and only learns from samples that are uniform on the hypercube. The authors in [12] deal with the problem of adaptive learning of a sparse polynomial where the algorithm has the power to query a specific sample. In learning theory, there is a lot of work [8, 17, 14, 5] on learning real valued sparse polynomials where the polynomial is evaluated on boolean samples from a product distribution. The closest to our work is [17] that gives a sub quadratic algorithm for learning any sparse quadratic polynomial of independent boolean variables with biases $1/2$ in sub-quadratic time. In our case, the setup involves a non-linearity and arbitrary biases of boolean variables although bounded away from $1/2$.

2 Problem Statement

In this paper we consider the sparse quadratic logistic regression setting as in (1). The main question we are interested in is:

Problem 1**.**

*Given samples $\langle y_{i},\mathbf{x}_{i}\rangle$ distributed according to the logistic regression model in 1, when $f$ is a quadratic polynomial containing at most $s$ parities, then is it possible to learn the function $f$ in **sub-quadratic time $\mathbf{\tilde{o}(p^{2})}$ *? If possible, is there an algorithm that runs in linear time $\tilde{O}(p)$?

$\tilde{o},\tilde{O}$ are equivalent the usual $o,O$ notation except in that it could include terms that are independent of the dimension $p$ (like sparsity, number of samples etc.). In this paper, we consider this question in the simplest case, when the $X_{i}$’s are independent. For the boolean case where $\mathcal{X}\in\{+1,1\}$ we assume $\Pr[X_{i}=1]=p_{i}$. Otherwise, we assume that $X_{i}$ is distributed according to a discrete pmf on the support $\mathcal{X}$.

Since there are $\Theta(p^{2})$ possible choices of non-zero coefficients, a log-likelihood optimization algorithm requires runtime that is polynomial in $p^{2}$ for estimating all the coefficients.

A correlation test with a single variable $X_{i}$ checks if $\lvert\mathbb{E}[Y(X_{i}-\mu_{i})]\rvert>\epsilon$ or not by using the empirical mean to calculate the correlation value. In this work, we discover a surprising property of this simple test in relation to logistic regression on binary variables: We first note that if the logistic function is not present and when the function is observed directly, i.e., $Y=f(\mathbf{X})$, then the weak-support cannot be detected from correlation of $Y$ with individual binary variables. We discuss this in the next section and show how the non-linearity $\sigma$ helps in identifying relevant variables even in cases when it is not possible to identify using correlation tests.

3 Simple Example: Non-linearity helps in the Boolean case

Let $\mathcal{X}\in\{+1,-1\}$. We now show that having the both the non-linearity $\sigma(\cdot)$, and biases of the input binary variables that is not either 0,1 or 1/2 exactly, is crucial to correlation tests (and hence also our algorithm) to work. We demonstrate this through a simple example and follow it up with an informal discussion on why it happens. Suppose we consider the problem of identifying the relevant binary variables $x_{i},x_{j}$ out of $p$ dependent variables given the labeled samples $(Y,\mathbf{x})$, where $Y=f(\mathbf{x})$. Let $f$ be a simple quadratic function $C(x_{i}-\mu_{i})(x_{j}-\mu_{j})$. Clearly, all correlation tests of the form $\mathbb{E}\left[Y(x_{k}-\mu_{k})\right]$, where $\mu_{k}=\mathbb{E}[x_{k}]$ are identically [math] for all $k\in[p]$. Hence, correlation tests cannot be used to identify $x_{i}$ or $x_{j}$. We are not aware of an algorithm to identify the weak support in linear time for these classes of functions.

Now, consider the same function $f$ with the logistic regression setup in $(\ref{def:Y})$. Let $\gamma=1$. It can be shown that $\mathbb{E}\left[Y(x_{k}-\mu_{k})\right]=\mathbb{E}\left[\sigma(f)\left(x_{k}-\mu_{k}\right)\right]\neq 0$ for $k={i,j}$ and [math] otherwise, when $p_{i},p_{j}$ are away from $0,1$ and $1/2$. Except in the non-trivial case when the variables are close to being uniform Bernoulli random variables, the correlation test works. Let $x_{1},x_{2}$ be the two variables in $f(\mathbf{x})$. We plot $\mathbb{E}[\sigma(f)(X_{1}-\mu_{1})]$ for $C=20,\gamma=1$ for all $p_{1},p_{2}\in[0,1]$ in Fig. 1. It is clear that, when $p_{1},p_{2}\neq\{0,1,1/2\}$, the correlation is non-zero. If the distance from these points is larger, then the correlation is also higher.

The quadratic polynomial $f$ is passed through a non-linear function $\sigma(\cdot)$ and we get to observe a noisy version of it (a coin toss with that probability). Here, the non-linearity $\sigma$ helps identify the variables using a simple correlation test. If $\sigma(\cdot)$ is linear, as in the above example, it is not possible to detect using a correlation test. This is the key message. We briefly explore why non-linearity helps and why it would not work with uniform random variables. In the next section, we show that for any set of bias probabilities away from $0,1$ and $1/2$, the correlation test succeeds when the unique sign property holds for the polynomial $f$ almost surely for scaling parameter $\gamma\in[0,1]$.

3.1 Key Idea behind our results: Why does non-linearity help?

Suppose $\sigma(x)$ is the polynomial $x^{q}$ for some positive integer $q$. Let $v_{1}\ldots v_{m}$ be the distinct non-zero values taken by $f$. Consider the influence of a relevant variable $x_{i}$ on the value $v_{i}$: $g(v_{i})=\mathrm{Pr}(\mathbf{x}:f(\mathbf{x})=v_{i}\lvert x_{i}=1)-\mathrm{Pr}(\mathbf{x}:f(\mathbf{x})=v_{i}\lvert x_{i}=-1)$. For an irrelevant variable, all influences are [math]. Then the correlation with $(x_{i}-\mu_{i})$ can be shown to be proportional to: $\gamma^{q}\sum v_{i}^{q}g(v_{i})$. If the test fails, then this specific linear combination of influences is [math]. Therefore, one can have non-zero influences and still fail the test.

Let us consider a non-linearity $\sigma(\cdot)$ that is analytic over the real line possessing an infinite degree taylor expansion. Now, the correlation with $(x_{i}-\mu_{i})$ is an infinite degree polynomial in $\gamma$, and for the test to fail for all $\gamma$, every coefficient has to be zero or the test fails for only countable values of $\gamma$. The coefficient of the $q$-th degree term is indeed $\sum v_{i}^{q}g(v_{i})$. This means that the influence variables $g(v_{i})$ are constrained by infinite equations, one for every degree $q$. For the test to fail for all $\gamma$, all these equations must be identically [math].

The unique sign property makes sure that $m$ of those equations are full rank and that there exists an influence variable that is not identically zero. Therefore, the test cannot fail for all $\gamma$. The presence of a strong non-linearity and the scaling parameter (along with the USP property) forces all influences to zero which is not possible. We show this formally in Theorem 5.1.

4 Algorithms

In this section, we outline our broad approach and give the main algorithms we propose. First, we have the following definition of weak and strong supports:

Definition 1**.**

The weak support of $f$, denoted by $S$, is defined to be the set of variables that $f$ depends on, i.e., $S=\{x_{i}:\exists j~{}\mathrm{s.t}~{}(i,j)\in Q~{}\mathrm{or}~{}i\in L\}$. Similarly, we define strong support of $f$ to be $\{Q,L\}$.

In order to achieve time complexity linear in $p$, we divide the objective of learning $f$ into three tasks:

Learning the weak support. 2. 2.

Learning the coefficients of the linear and quadratic terms formed by the weak support.

The last step can be implemented using a log-likelihood optimization in $O(s^{2})$ dimensions. In order to achieve sub-quadratic (in $p$) time complexity, it is sufficient to identify the weak support quickly. In fact, it is easy to see once the weak support is identified, the runtime of the other steps will not depend on $p$ but only on $s$. Hence, we focus on identifying the weak support first.

4.1 A Linear Correlation Test for Binary Variables

We first consider the model in (1) for binary variables. We propose Algorithm 1 for identifying the weak support which is based on a simple correlation test. Essentially, Algorithm 1 tests if $\lvert\mathbb{E}[Y(x_{k}-\mu_{k})]\rvert>\epsilon$ or not. The weak support algorithm runs in time $O(p)$ (Here, $O(\cdot)$ subsumes dimension independent quantities ).

Sample complexity: Algorithm 1 works when the number of samples $n$ is large enough, given that $\varepsilon$ is properly chosen. Given an $\varepsilon$ for which the population version succeeds, it can be shown that, if the number of samples is above $8c\log(p)\left(\frac{1}{\varepsilon}\right)^{2}$ the algorithm sssucceed with polynomially small error probability for any $c>1$. This is due to the fact that the random variable $Y(X_{i}-p_{i})$ lies in $[-2,2]$ and Hoeffding’s inequality can be used to derive strong concentration results. Existence of an $\varepsilon>0$ is shown for the sigmoidal non-linearity while we characterize an $\varepsilon$ for the case when $\sigma$ is piece-wise linear.

4.2 A Non-linear Correlation Test for Variables with Finite Support

In the previous sections, we saw intuitively that the non-linearity in the bias of the target variable leaves traces that can be detected through a simple correlation test. We will see in the subsequent sections that this statement is not unconditional: There are certain functions where the correlation test cannot be used to find the relevant variables. In this section, we propose a generalization of the correlation test, which we term as a non-linear correlation test. Non-linear correlation test applies a random non-linear transform on the features before finding their correlation with the target variable. Hence, the test can exploit random non-linearity injected in the test, rather than solely relying on the non-linearity due to nature. Notice that when the variable is binary $\{+1,-1\}$, any nonlinear test is equivalent to some linear test. Therefore variables of interest should be non-binary. Thus, we consider the model in (1) for real variables with finite support on the real line.

Through simulations, first we identify a specific function where the linear correlation test fails to detect one of the variables in the weak support. Later in Section 7.3, we empirically show that the nonlinear correlation test can be used to infer the existence of that variable in the weak support with reasonable number of samples.

4.2.1 Introducing Randomness through Hashing

We consider variables $x_{i}$ that take value from a finite set from the alphabet $\mathcal{X}\subset\mathbb{R}$. Let $f$ be a quadratic polynomial over $x_{i}$’s as in (1). We observe $Y\sim\mathrm{Ber}(\sigma(f(x)))$ where $\sigma$ is a non linearity. The nonlinear correlation test we propose is as follows: First pass $x_{i}$ through a random non-linearity $Z=g(x_{i})$ and perform the correlation test:

[TABLE]

We choose the non-linearity $g$ to be a randomly chosen hash function (from a hash family) each of which maps the values in $\mathcal{X}$ to the set of integers in $[-U,U]$. Note that, $\mathrm{stddev}(Z)$ is the standard deviation of the transformed random variable $Z=g(x_{i})$ for a fixed $g$ over the marginal distribution of $x_{i}$. Further, the $\mathbb{E}[\cdot]$ is over $Y,\{X_{i}\}$ and the randomness in choosing the hash $g(\cdot)$. We finally compare the obtained correlation values with some threshold and decide to include it in the weak support. We give a finite sample algorithmic version of this in Algorithm 2.

4.2.2 A simple example:

Let $X_{1},X_{2},X_{3}\in\{-2,-1,1,2\}$. Consider the function $f=\alpha(X_{1}-\mu_{1})(X_{2}-\mu_{2})+\beta(X_{3}-\mu_{3})+\gamma$. Let Y be the binary variable $Y\in\{-1,1\}$ were $\mathbb{P}(Y=1)=1/(1+\exp(-f))$. Using a simulation-aided search, we identify a set of coefficients $\alpha,\beta,\gamma$ and a set of marginal distributions for $X_{1},X_{2},X_{3}$ such that the correlation $E[(X_{1}-\mu_{1})Y]$ is close to zero. See Appendix B for the selected parameter values. Later, we pass $X_{1}$ through a non-linearity $Z=g(X_{1})$ and perform the correlation test in (2) via Algorithm 2. We choose the non-linearity to be a random hash function that maps the values $\{-2,-1,1,2\}$ uniformly randomly to integers in the range $[-1000,1000]$. We observe that even for 500 samples, in 99 out of 100 trials, $X_{1}$ has the nonlinear correlation higher than every irrelevant variable. For the number of samples larger than 1000, $X_{1}$ always has the highest nonlinear correlation.

Later, we show empirically that a large fraction of the variables in the weak support can be detected with the non-linear correlation test when the number of samples is large in the experimental section.

5 Analysis of Correlation Tests - Boolean Case

5.1 Analysis: Learning the Weak Support

In this section, we consider the model in (1) with boolean variables. We show the following smoothed analysis result for Algorithm 1: under mild assumptions on $f$ (needs to satisfy the unique sign property), weak support can be identified in linear time (in $p$) for almost all scaling parameters $\gamma\in[0,1]$ using a simple correlation test: Compare $E[Y(X_{i}-\mu_{i})]$ with [math]. The test succeeds when $E[Y(X_{i}-\mu_{i})]$ is [math] when $Y$ is independent of $x_{i}$ and not zero otherwise.

The function $f$ depends on $r$ binary variables. Therefore, it can take at most $2^{r}$ values. Let the distinct non zero values that the function takes be $v_{1},v_{2}\ldots v_{m}$ (where $m\leq 2^{r}$). Let $v_{m+1}=0$. This value is used if the function takes value [math].

5.1.1 Sufficient conditions

Without loss of generality, let us assume that the function $f$ depends on the first $r$ variables, i.e. $\{x_{1}\ldots x_{r}\}$. We will give conditions on $f$ and $\mu_{i}$’s such that the correlation with the first variable $\mathbb{E}[Y(x_{1}-\mu_{1})]\neq 0$.

Definition 2**.**

Define $G=(V,E)$ to be a graph on $r$ vertices associated with the function $f$ such that an edge $\{i,j\}\in E$ if $(i,j)\in Q$. Let $C$ be the connected component in the graph $G$ containing the variable $x_{1}$.

Definition 3**.**

Function $f$ is said to have a unique sign property with respect to a value $v_{i}$ if there is no $v_{j}\neq v_{i}$ such that $v_{j}=-v_{i}$ and when $f=v_{i}$, it fixes the sign of all parities $\{x_{i}x_{j}\}_{(i,j)\in Q},\{x_{j}\}_{j\in L}$.

Another way of stating the unique sign property is by saying $\lvert f\rvert=\lvert v_{i}\rvert$ fixes the sign of all parities $\{x_{i}x_{j}\}_{(i,j)\in Q},\{x_{j}\}_{j\in L}$.

We show that the idealized test succeeds when the bias probabilities are away from $0,1$ and $1/2$. There are two mutually exclusive cases in the main result. In one case, we need the function to satisfy the unique sign property for one non-zero value and in the other case we need it to satisfy the unique sign property for all its values.

Theorem 5.1**.**

Let $p_{i}\in(\delta,1/2-\delta)\cup(1/2+\delta,1-\delta),~{}\forall i$ for some $\delta>0$. Let $Y,\mathbf{x}$ follow the model given in (1). Let $\sigma(\cdot)$ be the sigmoidal function. Then, $E[Y(X_{1}-\mu_{1})]=0$ when $Y$ is independent of $X_{1}$. Also, $E[Y(X_{1}-\mu_{1})]\neq 0$ for all but finite values of $\gamma\in[0,1]$ when $Y$ is dependent on $X_{1}$ and when one of the following conditions are satisfied:

$f$* has a linear term in $x_{1}$ or has a linear term $x_{i}$ which is in the same component $C$ in $G$ as $x_{1}$ . $\exists v_{i}\neq 0:$ such that $\lvert f\rvert=\lvert v_{i}\rvert$ implies a unique sign for all the parities $\{x_{i}x_{j}\}_{(i,j)\in Q},\{x_{j}\}_{j\in L}$.* 2. 2.

$f$* has no linear term in $x_{1}$ and $C$ contains no variable with a linear term. $\forall$ values $v_{i}\neq 0:\lvert f\rvert=\lvert v_{i}\rvert$ implies a unique sign for all the parities $\{x_{i}x_{j}\}_{(i,j)\in Q},\{x_{j}\}_{j\in L}$.*

Proof.

The proof is relegated to the Appendix. The proof follows the basic outline provided in Section 3.1. ∎

Biases need to be away from $1/2$: We point out that if all bias probabilities are $1/2$, under the conditions in Case $2$, the test does not succeed. This is because the quantity that appears in the proof, i.e. $g(v_{i})$ is [math] for all $v_{i}$. This means that the bias probabilities being away from $1/2$ is essential for the final result.

Note: The condition that all the coefficients are in general position, i.e. $cs_{c}+\sum\limits_{(i,j)\in Q}\beta_{i,j}s_{i,j}+\sum\limits_{j\in L}\alpha_{j}s_{j}\neq 0$ for any set of signs $s_{i,j},s_{j},s_{c}\in\{0,1,-1\}$, ensures the unique sign condition stated for all values of the function in the Case $2$. For Case $1$, we have used the fact that there is no $j:v_{j}=-v_{i}$. This can be ensured if the constant $c$ is in general position with respect to other coefficients $\{\beta_{i,j}\},\{\alpha_{j}\}$. In fact, we only need, $c\neq\sum\limits_{(i,j)\in Q}\beta_{i,j}s_{i,j}+\sum\limits_{j\in L}\alpha_{j}s_{j}$ for any set of signs $s_{i,j},s_{j}\in\{0,1,-1\}$.

Remark: The above theorem works for a broad class of non-linearities $\sigma(\cdot)$ that is complex analytic and has more than $2^{s}$ non zeros terms in its taylor expansion. Any transcendental function is a special case that satisfies our properties.

6 Weak Support Learning: Finite Sample Complexity Results

In this section, we again consider boolean variables. We will show that when $Y$ depends on $x_{1}$ and when $f$ satisfies some technical conditions , $\lvert\mathbb{E}[Y(x_{1}-\mu_{1})]\rvert\geq\varepsilon$ for values of $\gamma$ in a set of measure which is significant in an interval on the real line. Here, we characterize $\varepsilon$ as a function of other parameters. This establishes sample complexity requirements for Algorithm 1. We prove this result when the non-linearity $\sigma(\cdot)$ is a piecewise linear function.

[TABLE]

The rest of the assumptions remain the same as in (1). Let $v_{1},\ldots v_{m}$ be the non zero values the function $f$ takes. For the purpose of the following theorem, let $\delta^{{}^{\prime\prime}}=\delta^{2r}\min(2^{\delta^{\prime}/2}-1,1-2^{-\delta^{\prime}/2})$ where $\delta^{\prime}=\left\lvert\log\left(\frac{2}{1-2\delta}\right)\right\rvert$ for some constant $\delta>0$.

Theorem 6.1**.**

Let the non-linearity $\sigma$ be the function in $(\ref{nonlinear})$. Assume that the coefficients of the function $f$ satisfy the following stronger general position property: $\lvert\sum\limits_{(i,j)\in Q}\beta_{i,j}s_{i,j}+\sum\limits_{j\in L}\alpha_{j}s_{j}+cs_{c}\rvert>\epsilon$ for all set of signs $s_{i,j},s_{j},s_{c}\in\{0,+1,-1\}$. Further, let the bias probabilities $p_{i}\in(\delta,1/2-\delta)\cup(1/2+\delta,1-\delta)$. When the function is independent of $x_{1}$, $\lvert\mathbb{E}[Y(X_{1}-\mu_{1})]\rvert=0$. When the function is dependent on $x_{1}$, we have the following cases:

If the component $C$ containing $x_{1}$ in graph $G$ has a linear term in $f$, then $\lvert\mathbb{E}[Y(X_{1}-\mu_{1})]\rvert>C_{1}\frac{\epsilon^{2}\delta^{r+2}}{b^{2}s^{2}}$ for $\gamma$ in a set of measure $C_{2}\frac{m\epsilon}{b^{2}s^{2}}$ in the open interval $(0,\frac{1}{\lvert v_{m}\rvert})\subseteq(0,\frac{1}{\epsilon})$, where $C_{1}=\frac{1}{32}$ and $C_{2}=\frac{3}{8}$ 2. 2.

If the component $C$ containing $x_{1}$ in graph $G$ has no linear term in $f$, then $\lvert\mathbb{E}[Y(X_{1}-\mu_{1})]\rvert>C_{1}\frac{\epsilon^{2}\delta^{2}\delta^{{}^{\prime\prime}}}{b^{2}s^{2}}$ for $\gamma$ in a set of measure $C_{2}\frac{m\epsilon}{b^{2}s^{2}}$ in the open interval $(0,\frac{1}{\lvert v_{m}\rvert})\subseteq(0,\frac{1}{\epsilon})$, where $C_{1}=\frac{1}{32}$ and $C_{2}=\frac{3}{32}$.

Proof.

A full proof appears in the Appendix. We give a brief sketch of the proof. The correlation with a relevant variable $x_{i}$ results in a piecewise linear function in the scaling parameter $\gamma$. Therefore, if the slope of most of these piecewise linear parts can be shown to be large, then a significant portion of the curve is away from the axis and implies a lower bound on the absolute value of the correlation. Although $\sigma$ is not analytic, we get a Vandermonde-like system of equations connecting influence variables ($g(v_{i})$ as in Section 3.1) to the slopes of the piecewise linear parts. Using a looser version of the unique sign property, one can show that a significant fraction of the influence variables are large in magnitude proving the result. ∎

7 Experiments

7.1 Synthetic Experiments: Weak Support Recovery

We illustrate the effectiveness of the correlation test in Algorithm 1 using synthetic experiments run on models that obey the setup in (1). We describe the experimental setup in Fig. 2. We fix the quadratic polynomial $f$ where $\lvert L\rvert$ linear terms are uniformly chosen from a fixed set of $r$ variables. Similarly, $\lvert Q\rvert$ bilinear terms are chosen uniformly randomly from the terms formed by $r$ variables. The coefficients of the quadratic polynomial is chosen from a uniform distribution over $[0.1,1]$. The number of relevant variables is $r=20$ and $s\sim 15-25$ for the experiments. For this fixed polynomial, we randomly generate $\gamma$ uniformly sampled from $[1,15]$ and for every gamma we generate several bias-probabilities uniformly randomly distributed in $[0.1,0.4]\bigcup[0.6,0.9]$. Here, $\delta$ is taken to be $0.1$. For a fixed polynomial, we average over $~{}140$ random initializations of $\gamma$ and bias-probabilities. We average over $20$ iterations with random choices of polynomials. We compute the fraction of times the weak support is recovered which is the ratio of the size of the part of the support set recovered to the the total size of weak support

Key Observations: The $90\%$ line for $R=20$ variables refers to the event: All but at most $2$ variables in the weak-support are recovered. In Fig. 2, we find that for different ratios of $\lvert L\rvert$ and $\lvert Q\rvert$, roughly $~{}10^{4}$ samples are required to cross this threshold when $p$ is in a few thousands. This is much smaller than the finite sample complexity results in Theorem 6.1. This suggests that in practice, the correlation tests work with fewer samples than what is predicted by the theory.

7.2 Experiments with real data-sets: Weak Support Recovery

Objective: We evaluate the effect of adding quadratic terms formed by the binary features in the weak support (of small fixed sparsity) identified by Algorithm 1 to existing binary features in a real data set. We show that this boosts the performance of a standard classifier that works with a linear function of the binary features (standard logistic regression and Linear SVM classifier).

Dataset: We performed our experiments on the Dorothea data set obtained from the UCI Machine Learning Repository. 111The URL for this data set is: https://archive.ics.uci.edu/ml/datasets/Dorothea. This data set was contributed to this repository by DuPont Pharmaceuticals Research Laboratories and KDD Cup 2001. This data set arises from a use case in drug design. In drug design, one wants a compound that can efficiently bind to a target receptor for the drug to be effective. The data set contains 1150 samples (800 train + 350 validation samples). Each sample is a set of $100000$ binary features of a compound. The data set has about $50000$ random features added to $50000$ genuine ones. The target is a binary variable indicating if the compound binds or not. We used the training samples for cross validation for our models while we tested on the validation set and we report the test AUC (Area under the Receiver Operating Characteristic Curve (ROC)) scores.

Algorithms Compared: We first perform a correlation test using Algorithm 1 and then we add all possible quadratic terms arising from the output weak support222When we do a correlation test using Algorithm 1, we normalize the correlation value by the standard deviation of that binary feature. Then, we rank the obtained normalized correlation values and pick the top $k$ features. to the existing features. We train a standard classifier that works with these expanded set of features. We compare the performance of the same classifier run without any feature addition. The regularization parameters $C$ in all our experiments used for cross validation is such that $\log_{10}(C)$ takes $15$ uniformly space values in the interval $[-4,4]$. We do a $4$ fold cross-validation at the training stage. More specifically, we compare the following four algorithms:

LSVM + $\ell_{1}$ regularization: This is a standard linear support vector machine with $\ell_{1}$ regularization. 2. 2.

LSVM + $\ell_{1}$ regularization + quadratic terms: Quadratic terms arising out of the $k$ top features according to the correlation test are added. For cross-validation, $k$ is chosen to lie in $\{3,5,8,10,15\}$. 3. 3.

*Logistic Regression + $\ell_{1}$ regularization * This is the standard logistic Regression with $\ell_{1}$ regularization. 4. 4.

Logistic Regression + $\ell_{1}$ regularization + quadratic terms: Quadratic terms arising out of the $k$ top features according to the correlation test are added. For cross-validation, $k$ is chosen to lie in $\{3,5,8,10,15\}$.

Results: We summarize the results in Table 1. We find a $0.02$ improvement in the AUC score for logistic regression and a $0.06$ improvement in the AUC score for Linear SVM. The running times with or without feature addition are comparable. We like to note that adding extra features to an already large feature matrix is a time consuming operation and time required for this step is included in the timing analysis.

7.3 Synthetic Experiments: Nonlinear Correlation Test for Variables with Finite Support

In this section, we provide empirical evidence that the non-linear correlation (Algorithm 2) identifies the weak support of a function dependent on real variables with finite support. We consider the following randomly chosen function: $f=\sum_{(i,j\in S\times S)}\alpha_{i,j}X_{i}X_{j}$, where $X_{i}\in\mathcal{X}=\{-2,-1,1,2\},\forall i$. The uniformly randomly chosen subset $S\subseteq[p]$ is the weak support of $f$, and $\alpha_{i,j}$ are i.i.d samples from a uniform random variable in the interval $[-1,1]$. The probability mass function of each variable $X_{i}$ is chosen uniformly randomly over the simplex in 4 dimensions. Hence, $f$ has all the quadratic terms, where the coefficient of each term is selected uniformly and independently in the range $[-1,1]$. We consider the case $p=1010$, and $\lvert S\rvert=10$. We perform the hashing-based nonlinear correlation test to obtain $C_{i}$ given in Section 4.2.1. Later, we choose the candidate weak support set $T$ by including every variable $X_{i}$ which is among the top 20 variables with highest correlation values $C_{i}$. We compute the fraction of the true weak support that is contained in $T$, i.e., $\frac{\lvert S\cap T\rvert}{\lvert S\rvert}$. This fraction is averaged over 100 randomly sampled functions and datasets, and is reported in Table 2 for varying number of samples. As observed, with enough samples, we can recover a large fraction of the variables in the weak support on average. Increasing the number of used hash functions improves the performance, although it incurs some computational cost.

8 Conclusion

We propose correlation tests to recover the set of relevant variables for the sparse quadratic logistic regression problem with real covariates with finite support. When the variables are all binary, the correlation test is a simple linear correlation test. We show that the non-linearity inherent in the problem helps the correlation test to succeed. Further, we propose a nonlinear correlation test that involves transforming covariates through hashing before performing correlation for the non-binary case. We show the effectiveness of our methods through a number of theoretical as well as experimental results.

Appendix A Appendix

A.1 Proof of Theorem 5.1

This section provides the proof of the theorem. We first develop some notation, state some intermediate Lemmas and then proceed to the main part of the proof.

It can be easily seen that:

[TABLE]

Let $E_{0}=1$ if the function takes value [math]. Let $\mathrm{Pr}(\mathbf{x}:f(\mathbf{x})=v_{i}|x_{k}=+1)-\mathrm{Pr}(\mathbf{x}:f(\mathbf{x})=v_{i}|x_{k}=-1)=g(v_{i}))$.

[TABLE]

The last equality is because $p_{k}(1-\mu_{k})=(1+\mu_{k})(1-p_{k})=2p_{k}(1-p_{k})$.

Let us assume that the non-linearity is sigmoidal as in $(\ref{def:Y})$. This means that $\sigma(x)=\frac{1}{2}+\frac{1}{2}\tanh\left(x/2\right)$. There is a domain $U$ (an open set) containing the real line in the complex plane such that $\tanh(\cdot)$ is an analytic function over $U$ [7]. This implies that $\left[\sigma(0)g(0)\mathbf{1}_{E_{0}=1}+\sum\limits_{i=1}^{m}\sigma(\gamma v_{i})g(v_{i})\right]$ is an analytic function in a domain $U^{\prime}$ (finite linear combinations of scaled version of analytics functions is analytic in an open neighborhood) containing the real line in the complex plane. It is well known that:

Theorem A.1**.**

[7]** If D is a domain on the complex plane and f(z) is a complex analytic function on $D$ that is not identically zero, then the zeros of $f$ in $D$ are isolated.

Theorem A.2**.**

[3]** (Littlewood-Offord Theorem) Consider the linear sum $\sum\limits_{k\in[n]}\alpha_{i}b_{i}$ where the joint distribution of $\{b_{i}\}$ is uniform on $\{-1,1\}^{n}$ and $\alpha_{i}\in\mathbb{R}$. If $\lvert\alpha_{i}\rvert>\nu$, then

[TABLE]

A set $E$ is said to be isolated if for every point $p\in E$, there is an open ball around $p$ that contains no element from $E$ other than p. This implies on the real line when $z=\gamma\in\left[0,1\right]$, the following is true:

Lemma A.3**.**

When $\sigma(\cdot)$ is an analytic in a domain containing the real line, exactly one of the following is true:

$\left[\sigma(0)g(0)\mathbf{1}_{E_{0}=1}+\sum\limits_{i=1}^{m}\sigma(\gamma v_{i})g(v_{i})\right]=0,~{}\forall\gamma\in[0,1]$. 2. 2.

$\left[\sigma(0)g(0)\mathbf{1}_{E_{0}=1}+\sum\limits_{i=1}^{m}\sigma(\gamma v_{i})g(v_{i})\right]=0$* for finitely many $\gamma$’s in $[0,1]$.*

Here, all the variables ($\gamma$ and $v_{i}$) and functions ($g$) are defined in $(\ref{eq:sigmoidsum})$.

**Note:**If the first condition is avoided, then clearly, $\mathbb{E}[Y(x_{k}-\mu_{k})]\neq 0$ a.s for all $\gamma\in[0,1]$. Clearly, when $Y$ is independent of $x_{k}$, then the first condition is true.

We now derive some sufficient conditions on $f$ and the bias probabilities $\{\mu_{i}\}$ such that the first condition is not true when $Y$ depends on $x_{k}$. This essentially proves that the idealized test works.

Proof of Theorem 5.1.

From Lemma A.3, it is enough to make sure that the first condition is not true. Since $\sigma(\cdot)$ is analytic in a domain containing the real line, there is a small interval $(0,\gamma_{0})\in[0,1]$ (with $\gamma_{0}\leq 1$) such that $\sigma(\cdot)$ has a taylor expansion around [math]. Keeping the sigmoidal function (an odd function with a constant shift) in mind, let us for simplicity assume that the taylor expansion is given by: $\sigma(x)=d_{0}+\sum\limits_{i=0}^{\infty}d_{2i+1}x^{2i+1}$ where $d_{0}\in\mathbb{R}\backslash 0$ and $d_{2i+1}\in\mathbb{R}\backslash 0,~{}\forall i$. Therefore, in the open interval $(0,\gamma_{1}=\min(\min\limits_{v_{i}\neq 0}\frac{\gamma_{0}}{\lvert v_{i}\rvert},\gamma_{0}))$, we have the following expansion:

[TABLE]

Suppose, the first condition in Lemma A.3 is true. Then clearly, $\sum\limits_{i=1}^{m}\sigma(\gamma v_{i})g(v_{i})=0$ in the interval $(0,\gamma_{1})$. Therefore, all the terms in $(\ref{taylorexp})$ are zero. $\sum\limits_{i=1}^{m}g(v_{i})+g(0)\mathbf{1}_{\{f\mathrm{~{}takes~{}}0\}}=\mathbb{E}[x_{1}-\mu_{1}]=0$ always. Therefore, the hypothesis that the first condition in Lemma A.3 is true implies:

[TABLE]

All the $v_{i}$’s are distinct. However, there could be $i\neq j:v_{i}=-v_{j}$. It can be easily seen that the number of linearly independent columns is $m-\lvert\{(i,j):v_{i}=-v_{j}\}\rvert$ for any finite truncation with $>m$ row. This is due to two reasons: 1) The matrix truncated is equivalent to a Vandermonde matrix when each column is diivded by some constant, i.e. $i$th column divided by $v_{i}$ 2) Since, $v_{i}$’s are distinct in their values, there is exactly one pair involving $i$ where $v_{i}=-v_{j}$ and only this causes reduction in rank.

For every pair $(i,j):v_{i}=-v_{j}$, because of the odd-power progression, the null space contains the following the vector with two ones: $[0~{}0~{}0~{}\underset{\mathrm{pos.}~{}i}{1}~{}\cdots~{}\underset{\mathrm{pos.}~{}j}{1}~{}\cdots~{}0]^{T}$. All these null-space vectors have disjoint support and hence orthogonal and linearly independent. This means that the null space is $\mathrm{span}\{[0~{}0~{}0~{}\underset{\mathrm{pos.}~{}i}{1}~{}\cdots~{}\underset{\mathrm{pos.}~{}j}{1}~{}\cdots~{}0]^{T}:\exists(i,j)~{}\mathrm{with~{}}v_{i}=-v_{j}\}$. Hence, the vector $[g(v_{1})\ldots g(v_{m})]^{T}$ avoid the null space of the matrix in $(\ref{vandeqn})$ when

[TABLE]

Under the condition in $(\ref{suffcond1})$, according to Lemma A.3, except for finite number of values of $\gamma\in[0,1]$, the idealized test succeeds since $[g(v_{1})\ldots g(v_{m})]^{T}$ avoids the null-space of the matrix in (20). Now, we give some conditions on the function $f$ and the probabilities $p_{i}$ of the variables when the above condition is satisfied. We will assume that $p_{i}\in(\delta,1/2-\delta)\cup(1/2+\delta,1-\delta),~{}\forall i$ for some $\delta>0$ in the all the cases that follow.

Case 1: $f$ has a linear term in $x_{1}$. $\exists v_{i}\neq 0:$ such that $\lvert f\rvert=\lvert v_{i}\rvert$ implies a unique sign for all the parities $\{x_{i}x_{j}\}_{(i,j)\in Q},\{x_{j}\}_{j\in L}$.

The signs of all the parities, uniquely determine the function value. Therefore, the above condition implies that there is no $v_{j}:v_{j}=-v_{i}$. From condition $(\ref{suffcond1})$, we just need to check if $g(v_{i})\neq 0$. Since every parity is fixed when $\lvert f\rvert=\lvert v_{i}\rvert$, , then $x_{1}$ is also fixed because $x_{1}$ has a linear term in $f$. Without loss of generality let us assume that when $f=v_{i}$, the sign of $x_{1}$ is fixed to be $1$. Based on (8) this implies, $g(v_{i})=\mathrm{Pr}(\mathbf{x}:f(\mathbf{x}=v_{i})|x_{1}=+1)$. Since all the bias probabilities $p_{i}$ are in $(\delta,1-\delta)$, $g(v_{i})\neq 0$.

Case 1’: $f$ has only quadratic terms in $x_{1}$ but there exists a linear term $x_{i}$ which is in the same component $C$ as $x_{1}$ in graph G. $\exists v_{i}\neq 0:$ such that $\lvert f\rvert=\lvert v_{i}\rvert$ implies a unique sign for all the parities $\{x_{i}x_{j}\}_{(i,j)\in Q},\{x_{j}\}_{j\in L}$.

Since $x_{i}$ and $x_{1}$ lie in the same component $C$ in graph $G$, there is a path $P\subseteq E$ from $x_{i}$ to $x_{1}$ containing adjacent edges in sequence. Since every parity is fixed when $f=v_{i}$, , then $x_{i}$ is also fixed because $x_{i}$ has a linear term in $f$. Let $x_{i},x_{j}$ be the first edge on $P$ from $x_{i}$. Since $x_{i}x_{j}$ is fixed and $x_{i}$ is also fixed, $x_{j}$ is also fixed. By propagating the signs across the path to $x_{1}$, the sign of $x_{1}$ is also fixed. Without loss of generality, let $x_{1}=+1$ when $f=v_{i}$. By the argument for the previous case, $g(v_{i})=\mathrm{Pr}(\mathbf{x}:f(\mathbf{x}=v_{i})|x_{1}=+1)$ and since all bias probabilities $p_{i}$ are $\delta$ away from [math] and $1$, $g(v_{i})\neq 0$.

Case 2: $f$ has only quadratic terms in $x_{1}$ and the connected component $C$ containing $x_{1}$ does not have any linear term in $f$. $\forall$ values $v_{i}\neq 0:\lvert f\rvert=\lvert v_{i}\rvert$ implies a unique sign for all the parities $\{x_{i}x_{j}\}_{(i,j)\in Q},\{x_{j}\}_{j\in L}$.

The assumption on the function leads to two properties: 1) For all $v_{i}$, $f=v_{i}$ implies a unique sign of the parities corresponding to $Q$ and $L$. 2) For all $v_{i}$, $\nexists v_{j}:v_{j}\neq-v_{i}$. Since the component of $x_{1}$ contains no linear terms, $x_{1}$ being $+1$ or $-1$ could result in the same sign for all parities. Because of the second property, we just need to show that $g(v_{i})\neq 0$ for some $i$. Consider a specific $v_{i}$. Due to the first property, the signs of parities in $Q$ and $L$ are fixed. Let value of the parity $x_{i}x_{j}$ with $(i,j)\in Q$ be $w_{i,j}$ while the value of $x_{j}$ with $j\in L$ be $y_{j}$. Now, consider the graph $G$ and the connected component containing $x_{1}$. When $x_{1}=+1$ and when $f=v_{i}$, every variable in the component $C$ takes a specific sign. This is because all parities are fixed. Every path from $x_{1}$ to another variable $x_{i}$ fixes $x_{i}$. This can be seen using arguments identical to Case $2$. For $x_{\ell}\in C$, let the sign taken by $x_{\ell}$ when $x_{1}=1$ be $z^{+}_{\ell}$. When $x_{1}=-1$, all the variables $x_{i}$ in $C$ change sign. Let the signs when $x_{1}=-1$ be $z_{\ell}^{-}=-z_{\ell}^{+}$. The other, variables outside $C$ are unaffected by the sign change of $x_{1}$. For ease of notation, define $q_{\ell}^{+}=(z^{+1}_{\ell}+1)/2$ and $q_{\ell}^{-}=(z^{-}_{\ell}+1)/2$. This means that $q^{+}_{\ell}=0\Rightarrow q^{-}_{\ell}=1$ and vice-versa. Therefore,

[TABLE]

Let,

[TABLE]

Now, $\lvert e(v_{i})\rvert*\left(\prod\limits_{j\in C-\{1\}}p_{j}\right)\geq\delta^{r}$ since $(1-p_{i}),p_{i}\geq\delta,~{}\forall i$. Therefore,

[TABLE]

$\frac{1-p}{p}$ is a decreasing function in $p$. Consider $\log\left(\frac{1}{p}-1\right)$ (with base $e$). At $p=1/2,~{}\log(1/p-1)=0$. Further, since $p_{i},~{}\forall i$ is away from $0,1$ and $1/2$ by $\delta>0$, $\infty>\log(1/p_{i}-1)>\log(\frac{2}{1-2\delta}-1)\mathrm{~{}or~{}}-\infty<\log(1/p_{i}-1)<\log(\frac{2}{1+2\delta}-1)$. Let $\delta^{\prime}=\lvert\log(\frac{1+2\delta}{1-2\delta})\rvert$. Then,

[TABLE]

Therefore, $(\ref{gbound1})$ becomes:

[TABLE]

Now, the variables $\left(1-2q^{+}_{j}\right)$ are binary variables which take values in $\{1,-1\}$. $v_{i}$ decides the parities uniquely and hence $q^{+}_{j}$ is unique. $f$ takes values $v_{i}$ and possibly [math]. There are totally $n=m$ or $n=m+1$ values depending on whether $f$ takes the value [math] or not. Suppose a value is drawn among the $n$ values randomly uniformly, we show that this induces a joint distribution of $\{q^{+}_{j}\}_{j\in C-\{1\}}$ which is uniform.

The component $C$ does not contain any linear term. Therefore, the signs of parities in $C$ and the rest of the parities are independent. Consider a spanning tree $T$ of $C$. It has $|C|-1$ (bilinear) parities and it has $|C|$ variables. When a parity $x_{i}x_{j}=+1$, then it corresponds to the linear equation $y_{i}+y_{j}=0$ over the binary field $\mathbb{F}_{2}$ where $y_{i}$ and $y_{j}$ are $0,1$ variables corresponding to $x_{i}$ and $x_{j}$ respectively. $+$ is an XOR operation over the field. $1$ is mapped to [math] and $-1$ is mapped to $1$. Now, the set of equations corresponding to the parities in $T$ ($|C|-1$ in number) are linearly independent over the binary field $\mathbb{F}_{2}$. Let the linear system be represented as $\mathbf{T}\mathbf{y}$ where $\mathbf{T}\in\mathbb{F}_{2}^{|C|-1\times|C|}$. Let us add an extra equation that sets $x_{1}$. This means that $y_{1}$ is set to a value. Now, the new set becomes $\mathbf{T}^{\prime}\mathbf{y}$ where the last row has a $1$ in the $1$st position corresponding to an equation for $y_{1}$. This corresponds to a full rank set of equations in $|C|$ variables. Now, if $\mathbf{T}^{\prime}\mathbf{y}=\mathbf{b}$, then $\mathbf{y}=(\mathbf{T}^{\prime})^{-1}\mathbf{b}$ because $\mathbf{T}^{\prime}$ is full rank. Now, if the parity values of $T$ along with the value of $x_{1}$ take uniform random values then the value of every variable in $C$ is also uniform. Clearly, when $x_{1}$ is fixed to be $1$, the remaining variables are also uniform random. This means that if the parity values in $C$ are uniform random then clearly, the $y_{i}$ values are uniformly random. This means that joint distribution of $\{q^{+}_{j}\}_{j\in C-\{1\}}$ is uniform.

Now, we show that uniform distribution on the $n$ values of $f$ induces a uniform distribution on the parities of $C$. Every value has a unique sign for the set of parities in $Q,L$. The set of parities in $C$ are independent of the rest of the parities. Therefore, the number of function values when the parities in $C$ are fixed is exactly $\frac{n}{2^{|C|-1}}$ and is exactly the same irrespective of the values of parities for $C$. This also means $n$ is a power of $2$. Therefore, a uniform distribution on the values will induce a uniform distribution on the parities of $C$ which has been shown to induce a uniform distribution on $\{q^{+}_{j}\}_{j\in C-\{1\}}$.

By $(\ref{lowerbndlog})$, Theorem A.2 and the arguments above with respect to a uniform random variable $V$ on the $n$ values of the function, we have:

[TABLE]

This means that for more than $\frac{1}{2}$ of the $v_{i}$’s (one less than this because $f$ could possible take [math]), $g(v_{i})\neq 0$.

∎

A.2 Proof of Theorem 6.1

Let the non zero function values be $\{v_{i}\}_{i=1}^{m}$. Let $v_{m+1}=0$ to be used when the function takes value [math]. Let the event that the function takes [math] be called $E_{0}$ (a $0,1$) boolean variable. The strong general position property implies that $\lvert v_{i}\rvert>\epsilon$ and $\lvert v_{i}-v_{j}\rvert>\epsilon$ and $\lvert v_{i}+v_{j}\rvert>\epsilon$. Together, this means that $\lvert\lvert v_{i}\rvert-\lvert v_{j}\rvert\rvert>\epsilon$. Further, any value (non-zero or zero) of the function determines a unique sign for the parities in $\{Q,L\}$. For any non-linearity $\sigma(\cdot)$, (8) holds. Therefore,

[TABLE]

Again, $\mathbb{E}[x_{1}-\mu_{1}]=0=p_{1}(1-p_{1})\left[\mathbf{1}_{E_{0}=1}g(0)+\sum\limits_{i=1}^{m}g(v_{i})\right]$. Let $g^{\prime}(v_{i})=\mathrm{sgn}(v_{i})g(v_{i})$. Therefore, $(\ref{correxpand})$ becomes:

[TABLE]

Since, $L$ is a piece wise linear function, the above expression is also piecewise linear with respect to $\gamma\in[0,\infty)$. Without loss of generality, let us assume that $v_{i}$ is ordered according to the the index $i$ such that $\lvert v_{i}\lvert<\lvert v_{j}\rvert,~{}\forall i>j)$. In the piecewise linear function in $(\ref{piecewiselinear})$, the function is linear between $\frac{1}{\lvert v_{i}\rvert}$ and $\frac{1}{\lvert v_{i+1}\rvert}$ $\forall i\geq 1$ and it is also linear between [math] and $\frac{1}{\lvert v_{1}\rvert}$ and linear in $(\frac{1}{\lvert v_{m}\rvert},\infty)$. Now, in every linear part, the function can cross [math] at most once. Hence, if the slope is not too small, the function is larger in absolute except for a very small interval around this zero. We first analyze the number of intervals among $(0,\frac{1}{\lvert v_{1}\rvert})$ , $(\frac{1}{\lvert v_{i}\rvert},\frac{1}{\lvert v_{i+1}\rvert}),~{}1\leq i\leq m-1$ and $(\frac{1}{\lvert v_{m}\rvert},\infty)$ where the magnitude of the slope is large. Let $b_{0},b_{1}\ldots b_{m-1},b_{m}$ be the slopes of these intervals. Then, we have the following linear system from the definition of $L(\cdot)$:

[TABLE]

First, we observe that since $\lvert\lvert v_{i}\rvert-\lvert v_{j}\rvert\rvert>\epsilon$ for all $i\neq j$, there is at most one $\lvert v_{i}\rvert$ such that $\lvert\lvert v_{i}\rvert-1\rvert<\epsilon/4$. This means that, for all $i$ except one $\lvert\lvert v_{i}\rvert-1\rvert>\epsilon/4$.

Case 1: When the function $f$ has a linear term in $x_{1}$ or if the component $C$ containing $x_{1}$ has a linear term, by the same arguments in Theorem 5.1, $g(v_{i})=\mathrm{Pr}\left(\mathbf{x}:f(\mathbf{x})=v_{i}|x_{1}=+1\right)$ or $g(v_{i})=-\mathrm{Pr}\left(\mathbf{x}:f(\mathbf{x})=v_{i}|x_{1}=-1\right)$ for all $i$. Here, every values $v_{i}$ implies a unique sign for the parities. Clearly, since $p_{i},1-p_{i}\geq\delta$, in this case: $\lvert g^{\prime}(v_{i})\rvert\geq\delta^{r},~{}\forall i$. Now, consider the matrix equation in $(\ref{eqn:slope})$. If, for a particular $i$, $\lvert b_{i}\rvert<\frac{\epsilon}{8}\delta^{r}$, then both $\lvert b_{i-1}\rvert>\frac{\epsilon}{8}\delta^{r}$ and $\lvert b_{i+1}\rvert>\frac{\epsilon}{8}\delta^{r}$. This is because consecutive $b_{i}$’s differ by $(\lvert v_{i}\rvert-1)g^{\prime}(v_{i})$ which is bounded below in absolute value. This means that at least $m/2$ slopes among $b_{0},\ldots b_{m-1}$ are large in magnitude. Each interval is of length at least $\lvert\frac{1}{|v_{i}|}-\frac{1}{|v_{i+1}|}\rvert\geq\frac{\epsilon}{b^{2}s^{2}}$. It is easy to verify that in an interval $(\frac{1}{\lvert v_{i}\rvert},\frac{1}{\lvert v_{i+1}\rvert})$ whose slope is greater than $\frac{\epsilon}{8}\delta^{r}$, except for a smaller interval of length $\frac{\epsilon}{4b^{2}s^{2}}$ contained in it, the correlation $\lvert\mathbb{E}[Y(X_{k}-\mu_{k})]\rvert>\frac{\epsilon^{2}\delta^{r+2}}{32b^{2}s^{2}}$ (since $(1-p_{1})p_{1}\geq\delta^{2}$). Therefore, the test succeeds for $\gamma$ in a set of measure $\frac{3m\epsilon}{8b^{2}s^{2}}$ in the open interval $(0,\frac{1}{\lvert v_{m}\rvert})\subseteq(0,\frac{1}{\epsilon})$.

Case 2: Consider the case when $f$ is such that the component containing $x_{1}$ has no linear term. Further, from the assumptions, $f$ is such that $|f|=\lvert v_{i}\rvert$ fixes the signs of all parities. Therefore, the value of $g(v_{i})$ is identical to Case $2$ in Theorem 5.1. From $(\ref{eqn:halfbig})$, more than $\left(1-\frac{O(1)}{\sqrt{|C|}}\right)$ fraction of the $\lvert g^{\prime}(v_{i})\rvert$’s are larger than $\delta^{{}^{\prime\prime}}=\delta^{2r}\min(2^{\delta^{\prime}/2}-1,1-2^{-\delta^{\prime}/2})$ where $\delta^{\prime}=\left\lvert\log\left(\frac{2}{1-2\delta}\right)\right\rvert$.

Let ${\cal I}=\{i:\lvert g^{\prime}(v_{i})\rvert>\delta^{{}^{\prime\prime}}\}$. From the discussion in the preceding paragraph, $\lvert{\cal I}\rvert\geq m(1-\frac{O(1)}{\sqrt{|C|}})$. Now, clearly, for every $i\in{\cal I}$ either $b_{i}$ or $b_{i+1}$ is large ($>\delta^{{}^{\prime\prime}}\frac{\epsilon}{8}$) in magnitude because $\lvert|v_{i}|-1\rvert\lvert g^{\prime}(v_{i})\rvert\geq\delta^{{}^{\prime\prime}}\frac{\epsilon}{4}$. Therefore, there can be at least $\frac{\lvert\cal I\rvert}{4}$ slopes among $b_{0},\cdots b_{m-1}$ that is large. Therefore, modifying the constants in the previous case with constants derived above, we have: For $\gamma$ in a set of measure $\frac{3m\epsilon}{16b^{2}s^{2}}(1-\frac{O(1)}{\sqrt{|C|}})$ in the open interval $(0,\frac{1}{\lvert v_{m}\rvert})\subseteq(0,\frac{1}{\epsilon})$, $\lvert\mathbb{E}[Y(X_{k}-\mu_{k})]\rvert>\frac{\epsilon^{2}\delta^{2}\delta^{{}^{\prime\prime}}}{32b^{2}s^{2}}$.

Appendix B Parameter Values for the Function in Section 4.2.1

Following are the selected parameters: $\alpha=2.008,\beta=2,\gamma=2$. Probability mass function for the variables are given by

[TABLE]

Appendix C Strong Support Identifying Algorithms

We assume that $f$ depends on $r$ variables $\{x_{1},x_{2},...,x_{r}\}$. We assume that the weak support learning algorithm has identified the weak support containing the $r$ relevant variables. The objective now is to determine whether $x_{i}$ appears as a linear term, quadratic term, or both in the function. We study Algorithm 3 for an important special class of functions.

We have the following proposition regarding the above algorithm for a very special class of functions.

Proposition 1**.**

Assume that each variable $x_{i}$ appears once in the function $f$, i.e., either as $\alpha_{i}x_{i}$ for some $i\in L$ or as $\beta_{i,j}x_{i}x_{j}$ for some $i,j\in Q$ but not both. Then, $\mathbb{E}[Y|X_{i}=1,X_{j}=1]=\mathbb{E}[Y|X_{i}=-1,X_{j}=-1]$ and $\mathbb{E}[Y|X_{i}=1,X_{j}=-1]=\mathbb{E}[Y|X_{i}=-1,X_{j}=+1]$ when $(i,j)\in Q$ and one of the checks is not true when $(i,j)\notin Q$.

Proof.

In order to prove the theorem, following simple lemmas are useful:

Lemma C.1**.**

For any $f_{1}(x),f_{2}(x)$, we have

[TABLE]

Proof.

[TABLE]

The proof follows from the fact that $\exp(t)>0,\forall t$. ∎

Lemma C.2**.**

Let $c\geq 0$ be a constant and $f(x)$ be a bounded function. Then we have

[TABLE]

with equality if and only if $c=0$.

Proof.

For random variable $x\in[m]$, and $p_{i}=\Pr(x=i)$,

[TABLE]

If $c=0$, equality is trivial. If $c>0$, then $\sigma(\alpha+c)>\sigma(\alpha)$ for any $\alpha$ since $\sigma$ is a monotonically increasing function; hence the strict inequality. ∎

It is sufficient to prove that $S$ contains all the quadratic terms and no linear terms.

Case 1: $x_{i}$ and $x_{j}$ appears as a quadratic term, i.e. $\beta_{i,j}x_{i}x_{j}$:

Since $x_{i},x_{j}$ does not appear in another term, by the assumption of the function, expectation does not change unless the value of $x_{i}x_{j}$ changes, $\mathbb{E}[Y|X_{i}=1,X_{j}=1]=\mathbb{E}[Y|X_{i}=-1,X_{j}=-1]$ and $\mathbb{E}[Y|X_{i}=1,X_{j}=-1]=\mathbb{E}[Y|X_{i}=-1,X_{j}=+1]$

Case 2: $x_{i}$ and $x_{j}$ appear as linear terms, $\alpha_{j}x_{j}$ and $\alpha_{i}x_{i}$:

Then, for the first pair of checks in the test, we have $\mathbb{E}[Y|X_{i}=1,X_{j}=1]=\mathbb{E}[\sigma(f_{0}+\alpha_{i}+\alpha_{j})]$ and $\mathbb{E}[\mathbb{E}[Y|X_{i}=-1,X_{j}=-1]=\mathbb{E}[\sigma(f_{0}-\alpha_{i}-\alpha_{j})]$. Defining $h_{0}=f_{0}-\alpha_{i}-\alpha_{j}$, we have $\mathbb{E}[\sigma(h_{0}+2(\alpha_{i}+\alpha_{j}))]=\mathbb{E}[\sigma(h_{0})]$. Lemma C.2 implies $\alpha_{i}=-\alpha_{j}$. Similarly, second check implies $\alpha_{i}=\alpha_{j}$, which is not possible. Hence, the check fails in this case.

Case 3: $x_{i}$ appears in a linear term $\alpha_{i}x_{i}$ and $x_{j}$ appears in a quadratic term, $\beta_{j,k}x_{j}x_{k}$:

Let us assume that both checks succeed and show a contradiction. Expand $\mathbb{E}[\sigma(f+\alpha_{i}X_{i}+\beta_{j,k}X_{j}X_{k})]$ as

[TABLE]

First check yields (denoting $p(X_{k}=1)$ by $p$)

[TABLE]

and the second check yields

[TABLE]

Let $h(\mathbf{X})=f(\mathbf{X})-\alpha_{i}-\beta_{j,k}$. Then the equations become

[TABLE]

and

[TABLE]

respectively. Rearranging yields

[TABLE]

and

[TABLE]

Renaming the difference of expectation terms, we have $Ap=B(1-p)$ and $A(1-p)=Bp$. This implies that, $A=B$ and $p=0.5$. But $p=0.5$ is not possible. This contradicts the assumption that the checks succeed.

Case 4: $x_{i}$ and $x_{j}$ both appear as distinct quadratic terms, $\beta_{i,k}x_{i}x_{k},\beta_{j,\ell}x_{j}x_{\ell}$:

Similar to case 3, we assume the checks succeed and show a contradiction. Expand $\mathbb{E}[\sigma(f+\mathbf{\beta_{i,k}X_{i}X_{k}}+\mathbf{\beta_{j,\ell}X_{j}X_{\ell}}]$ by conditioning on both $X_{k}$ and $X_{\ell}$

[TABLE]

First check yields (denoting $p(X_{k}=1)$ by $p$, $p(X_{\ell}=1)$ by $q$, and $\bar{p}=1-p,\bar{q}=1-q$)

[TABLE]

and the second check yields

[TABLE]

Let $h(\mathbf{X})=f(\mathbf{X})-\beta_{i,k}-\beta_{j,\ell}$. Then the equations become

[TABLE]

and

[TABLE]

respectively. Rearranging yields

[TABLE]

and

[TABLE]

Renaming the difference of expectation terms, we have $Apq+Bp\bar{q}=B\bar{p}q+A\bar{p}\bar{q}$ and $Bpq+Ap\bar{q}=A\bar{p}q+B\bar{p}\bar{q}$. Summing both equations yield $(A+B)p=(A+B)(1-p)$, which implies either $A=-B$ or $p=0.5$. Assume $A=-B$. Then we have,

[TABLE]

Since $A=\mathbb{E}[\sigma(h+2\beta_{i,k}+2\beta_{j,\ell})]-\mathbb{E}[\sigma(h)]$ cannot be zero unless $\beta_{i,k},\beta_{j,\ell}$ are zero. Hence either $p=0.5$ or $q=0.5$, which is not allowed by assumption.

∎