@Randel, sorry to keep you waiting for over 5 years, but I just came up with a cool little hack that apparently solves for this. It's based on iteratively updating a constant offset until convergence. The R code below simulates data from a logistic model with $logit(\Pr(Y= 1 | x_1, x_2)) = -2 + .5 x_1 + .5x_2$, and creates a situation where both regression parameters are equal.

set.seed(124)
N <- 500
sim_df <- data.frame(x1 = rnorm(N))
sim_df$x2 <- 1.3 + .5 + sim_df$x1 + rnorm(N)
cor(sim_df$x1, sim_df$x2) # to keep it interesting

sim_df$z <- with(sim_df,  -2 + .5 * x1 + .5 * x2)
sim_df$u <- runif(N)

sim_df$y <- with(sim_df, as.numeric(u < plogis(z)))
table(sim_df$y)

my_model <- glm(y ~ x1 + x2, data = sim_df, family = "binomial")
summary(my_model)

Even though the true values of the coefficients are the same, the estimated values are clearly different:

Call:
glm(formula = y ~ x1 + x2, family = "binomial", data = sim_df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.8887  -0.8124  -0.5640   0.9315   2.5558

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.8947     0.2382  -7.956 1.78e-15 ***
x1            0.3558     0.1557   2.284   0.0224 *
x2            0.5467     0.1109   4.928 8.30e-07 ***
---

Let's initialize based on the above model and iterate:

beta1 <- coef(my_model)[["x1"]]
prev_beta1 <- 0

while (abs(prev_beta1 -  beta1) > .00001) {
  my_model <- glm(y ~ x1 + offset(x2 * beta1), data = sim_df,
                  family = "binomial")
  prev_beta1 <- beta1
  beta1 <- coef(my_model)[["x1"]]
  cat("previous beta1: ", prev_beta1, "new beta1: ", beta1, "\n")
}

  -- This is output:
  previous beta1:  0.3557504 new beta1:  0.5203537
  previous beta1:  0.5203537 new beta1:  0.3778784
  previous beta1:  0.3778784 new beta1:  0.5007586
  previous beta1:  0.5007586 new beta1:  0.3944498
  previous beta1:  0.3944498 new beta1:  0.4861758
  previous beta1:  0.4861758 new beta1:  0.406849
  previous beta1:  0.406849 new beta1:  0.4753156

summary(my_model)

Call:
glm(formula = y ~ x1 + offset(x2 * beta1), family = "binomial",
    data = sim_df)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.8289  -0.8224  -0.5777   0.9677   2.4802

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.6990     0.1077  -15.78  < 2e-16 ***
x1            0.4435     0.1239    3.58 0.000343 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 541.75  on 499  degrees of freedom
Residual deviance: 527.63  on 498  degrees of freedom
AIC: 531.63

Number of Fisher Scoring iterations: 4

Note that the final model with the offset doesn't include $x_2$, but you can see below that it's in the model, with the same coefficient value as $x_1$:

# compare:
predict(my_model)[1]
# vs
coef(my_model)[1] +
  coef(my_model)[["x1"]] * sim_df[1, "x1"] +
  coef(my_model)[["x1"]] * sim_df[1, "x2"]