How to recover linear probability model for binary classification? In some binary classification problem I assume that the probability for a positive is exactly linearly dependent on the features $P(y=1|x_1,x_2,x_3,\ldots)=\beta_1x_1+\beta_2x_2+\beta_3x_3+\cdots$. Now for the training I only observe the realization 0 or 1 for $y$. Therefore the probability isn't directly observable, but would rather be some local density (which is hard to estimate for high dimensions).
Is there some type of regression that can determine the coefficients $\beta_i$?
I could try normal regression, but it doesn't seem obvious that any loss function will be able to converge to the same $\beta_i$. What type of regularization and loss function would be most appropriate?
 A: This is the standard Linear Probability Model in the context of Binary Response models. To obtain it, one assumes that there exist an underlying latent/unobservable variable $y^*$, which is linearly dependent on the regressors,
$$y^*=\beta_1x_1+\beta_2x_2+\beta_3x_3+\cdots + u$$
and whose indicator function is the observed $y$, indicating that the latent variable exceeds (or falls below) some threshold (usually zero). In order to arrive at the linear specification for the probability, we have to assume that the error term $u$ in the underlying regression is uniformly distributed, conditional on the regressors.
See this answer for the exact derivation. The model can of course be estimated by least-squares, but has various drawbacks, so frequently a "probit" or "logit" model is used (the first assuming that the underlying error is normally distributed, the second assuming that the underlying error term follows the logistic distribution). These models are usually estimated by maximum likelihood.
ADDENDUM
After some discussion with the OP and clarification in the comments, it is clear that the OP's framework is represented by the linear probability model where the "error term" is the random draws the OP performs from a $U(0,1)$. Specifically we can postulate the unobservable regression
$$y^*_i = -\sum_{j=1}^k\beta_jx_{ji} + u_i,\qquad u_i \sim U(0,1)$$
and the OP determines the values of the observable $y$ by
$$y_i = 1,\;\;  \text{if} \;\;u_i \leq \sum_{j=1}^k\beta_jx_{ji} \Rightarrow - \sum_{j=1}^k\beta_jx_{ji} +u_i \leq 0 \Rightarrow y^*_i \leq 0$$
So we have
$$P(y_i = 1 \mid \mathbf X) = P\left(- \sum_{j=1}^k\beta_jx_{ji} +u_i\leq 0\right) =P\left(u_i\leq \sum_{j=1}^k\beta_jx_{ji} \right) $$
But this is the expression for the cumulative distribution function of $u_i$, which is $F_U(z) = P\left(u_i\leq z\right) = z$ for a $U(0,1)$.
In this way we arrive at the equation
$$P(y_i = 1 \mid \mathbf X) = \sum_{j=1}^k\beta_jx_{ji}$$
which is the Linear Probability Model, on firm stochastic ground. Note that now the "error term" is no longer "present" -but it has determined the structure of the equation, through its stochastic properties.  
Estimating the $\beta$'s from this last equation by Least Squares does not require specifying some additional "error term" -it has now become more of a mathematical approximation endeavor, rather than a stochastic regression. Another way of course would be to estimate the model by maximum likelihood. Here the distributional specification on $y$ is "safe": it is by construction a Bernoulli random variable, whose conditional parameter $p$ is determined by $\sum_{j=1}^k\beta_jx_{ji}$. We can write the Bernouli probability mass function as 
$$ f_Y(y_i\mid \mathbf X) = p_i^{y_i}(1-p_i)^{1-y_i} = \left(\sum_{j=1}^k\beta_jx_{ji}\right)^{y_i}\left(1-\sum_{j=_1}^k\beta_jx_{ji}\right)^{1-y_i}$$
from which we can proceed to form the log-likelihood of the sample and perform maximum likelihood.  
Estimator properties and other related issues to the estimation are easily accessible.
