Best approach to find optimal solution to linear equation by group in R
You could use linear regression without intercept in order to accomplish this.

If I understand correctly, your implementation of linear regression in Update 3 had two issues: first, you fit the linear model with intercept, but then drop the intercept when applying the model. Second, there is a typo in the following line

linear_discount = 1/linear_discount-1


In the code snippet below I extended your analyses by two approaches: linear regression without intercept and linear regression with intercept. The former case should correspond exactly to your requirements, i.e. using a single relative discount rate per customer. The latter one corresponds to the case with relaxed assumptions, i.e. one relative discount rate and one absolute discount (independent on the price level).

Note that for the illustrative dataset you provided the optimal (in terms of mean squared error) discount rate would be zero.

df_ext &lt;- df %&gt;%
  # mean hack  
  mutate(Discount_indiv = (Old_Price-New_Price)/New_Price) %&gt;% # Identify optimal discount individually
  group_by(Customers) %&gt;%
  mutate(Discount_mean = mean(Discount_indiv),
         Old_Price_mean = New_Price + New_Price*Discount_mean) %&gt;% # Average individual discount to get approximate discount for customer

  # incorrect linear regression (fitted with intercept, then intercept dropped)  
  do({ co &lt;- coef(lm(Old_Price ~ New_Price, .))
  mutate(., Discount_wrong = co[2])
  }) %&gt;%
  ungroup %&gt;%
  mutate(Discount_wrong = 1 / Discount_wrong-1) %&gt;%
  mutate(Old_Price_wrong = New_Price+New_Price*Discount_wrong) %&gt;%

  # correct linear regression w/o intercept
  group_by(Customers) %&gt;%
  mutate(Discount_regr = coef(lm(Old_Price ~ New_Price -1))[1] - 1,
         Old_Price_regr1 = New_Price+New_Price*Discount_regr) %&gt;% 

  # correct linear regression with intercept (i.e. two discounts per customer)
  mutate(
    Discount_abs = coef(lm(Old_Price ~ New_Price))[1],
    Discount_rel = coef(lm(Old_Price ~ New_Price))[2] - 1,
    Old_Price_regr2 = New_Price + New_Price*Discount_rel + Discount_abs) %&gt;% 

  #calculate residuals
  ungroup() %&gt;%
  mutate(Resid_mean = Old_Price_mean - Old_Price,
         Resid_wrong = Old_Price_wrong - Old_Price,
         Resid_regr1 = Old_Price_regr1 - Old_Price,
         Resid_regr2 = Old_Price_regr2 - Old_Price)


#transform data for visualisation
df_gat &lt;- select(df_ext, matches("Customers|Price")) %&gt;% 
  gather(key="Approach", value="Old_Price", -Customers, -New_Price)

ggplot(df_gat, aes(x=New_Price, y=Old_Price, group=Approach, colour=Approach, shape=Approach)) + 
  geom_line() + geom_point(size=3) + facet_wrap(~Customers, ncol = 1)

select(df_ext, matches("Customers|Disc"))
select(df_ext, matches("Customers|Price"))

#calculate mean square error for all approaches
select(df_ext, matches("Customers|Res")) %&gt;% 
  group_by(Customers) %&gt;%
  summarise_all(~mean(.^2))


Results for all approaches: MSE by customer &amp; some visualizations

  Customers Resid_mean Resid_wrong Resid_regr1 Resid_regr2
  &lt;fct&gt;          &lt;dbl&gt;       &lt;dbl&gt;       &lt;dbl&gt;       &lt;dbl&gt;
1 A               1.71        3.62       1.          0.222
2 B               1.71       11.5        1.000       0.222