java - Python scipy.optimize.leastsq 到 Java org.apache.commons.math3.fitting.leastsquares

Question

我尝试模仿这个algorithm ，用 Python 开发，根据看到的 Wifi 站位置计算地理位置，本身基于此 idea .

该算法首先使用 Numpy 函数来计算观测到的纬度和经度的基本加权平均值。为了最大限度地减少可能的 Wifi 位置错误的影响，它还使用“ scipy.optimize.leastsq ”方法，以便以统计方式计算，并在可能的情况下计算更精确的位置。

我想在 Java Android 平台上实现相同的行为。

对于所有其他计算，我成功地依赖于 org.apache.commons.math3。因此，对于最小二乘问题，我在逻辑上尝试依赖 https://commons.apache.org/proper/commons-math/userguide/leastsquares.html .

如果我很好理解的话，我的问题是 Scipy 为我管理雅可比函数定义的复杂性，而我糟糕的数学技能不允许我正确定义最小二乘问题的模型。我基于这个example尝试了一些实验，这似乎接近我需要的，但结果不好，因为我不知道如何处理“雅可比”部分。

就像有人为 this 所做的那样帖子，有人可以为我做同样的事情并尝试以简单的方式解释它吗？

有关 Python 部分如何工作的更多详细信息:

使用的“scipy.optimize.leastsq”语句是:

(lat, lon), cov_x, info, mesg, ier = 
scipy.optimize.leastsq(func, initial, args=data, full_output=True)

其中数据为:纬度/经度/以毫秒为单位的年龄/信号强度，例如:data = numpy.array([(43.48932915, 1.66561772, 1000, -20), (43.48849093, 1.6648176, 2000, - 10), (43.48818612, 1.66615113, 3000, -50)])

Initial 是计算出的加权平均纬度/经度，在此示例中:initial = 43.48864654, 1.66550075

功能是

  def func(initial, data):
        return numpy.array([
            geographic distance((float(point[latitude]), float(point[longitude])), (initial[latitude], initial[longitude])).meters * min(math.sqrt(2000.0 / float(point[age])), 1.0) / math.pow(float(point[signal strength]), 2)

结果是:43.4885401095, 1.6648660983

我在 Java 中进行的实验，我替换了数据值并更改了“modelI”的计算方式。我简化了信号强度和年龄值。但这是事实，而且结果表明，这还不够。

double modelI = calculateVincentyDistance(o.getY(), o.getX(), center.getY(), center.getX())* Math.min(Math.sqrt(2000.0/1000.0), 1.0) / Math.pow(-10, 2);

我也要去尝试https://github.com/odinsbane/least-squares-in-java ，但我不确定是否正确使用它，因为我不掌握它的工作方式。

仅供引用，我使用文森蒂距离计算，例如可以用半正矢或欧几里得代替。

感谢您的帮助!

Answer 1

代码不容易移植，因为 SciPy 提供了更通用的 Least-squares minimization接口(interface)，而 Apache Commons Math 提供 curve fitting 。还有许多优化问题可以重新表述为曲线拟合。在你最小化的Python代码中

F(current_point) = Sum{ (distance(known_point[i], current_point) * weight[i])^2 } -> min

Java 曲线拟合问题有点不同:

F(current_point) = Sum{ (target_value[i] - model[i](current_point))^2  } -> min

因此，可以通过将所有 target_value 分配为 0 并使 model[i] 计算从 current_point 到known_point[i]。

在一般情况下，此类问题没有使用公式的精确解，而是使用一些数值优化方法。这里还有另一个区别:Java 实现明确要求您为优化器提供计算正在优化的函数的导数的方法。如果未提供 Dfun，Python 代码似乎会使用某种差异区分器。您可以在 Java 中手动或使用 FiniteDifferencesDifferentiator 执行类似的操作但对于简单的公式，使用 DerivativeStructure 显式编码可能更容易

static class PositionInfo {
    public final double latitude;
    public final double longitude;
    public final int ageMs;
    public final int strength;

    public PositionInfo(double latitude, double longitude, int ageMs, int strength) {
        this.latitude = latitude;
        this.longitude = longitude;
        this.ageMs = ageMs;
        this.strength = strength;
    }

    public double getWeight() {
        return Math.min(1.0, Math.sqrt(2000.0 / ageMs)) / (strength * strength);
    }
}


static DerivativeStructure getWeightedEuclideanDistance(double tgtLat, double tgtLong, PositionInfo knownPos) {
    DerivativeStructure varLat = new DerivativeStructure(2, 1, 0, tgtLat); // latitude is 0-th variable of 2 for derivatives up to 1
    DerivativeStructure varLong = new DerivativeStructure(2, 1, 1, tgtLong); // longitude is 1-st variable of 2 for derivatives up to 1
    DerivativeStructure latDif = varLat.subtract(knownPos.latitude);
    DerivativeStructure longDif = varLong.subtract(knownPos.longitude);
    DerivativeStructure latDif2 = latDif.pow(2);
    DerivativeStructure longDif2 = longDif.pow(2);
    DerivativeStructure dist2 = latDif2.add(longDif2);
    DerivativeStructure dist = dist2.sqrt();
    return dist.multiply(knownPos.getWeight());
}

// as in https://en.wikipedia.org/wiki/Haversine_formula
static DerivativeStructure getWeightedHaversineDistance(double tgtLat, double tgtLong, PositionInfo knownPos) {
    DerivativeStructure varLat = new DerivativeStructure(2, 1, 0, tgtLat);
    DerivativeStructure varLong = new DerivativeStructure(2, 1, 1, tgtLong);
    DerivativeStructure varLatRad = varLat.toRadians();
    DerivativeStructure varLongRad = varLong.toRadians();
    DerivativeStructure latDifRad2 = varLat.subtract(knownPos.latitude).toRadians().divide(2);
    DerivativeStructure longDifRad2 = varLong.subtract(knownPos.longitude).toRadians().divide(2);
    DerivativeStructure sinLat2 = latDifRad2.sin().pow(2);
    DerivativeStructure sinLong2 = longDifRad2.sin().pow(2);
    DerivativeStructure summand2 = varLatRad.cos().multiply(varLongRad.cos()).multiply(sinLong2);
    DerivativeStructure sum = sinLat2.add(summand2);
    DerivativeStructure dist = sum.sqrt().asin();
    return dist.multiply(knownPos.getWeight());
}

使用这样的准备，你可以做这样的事情:

public static void main(String[] args) {

    // latitude/longitude/age in milliseconds/signal strength
    final PositionInfo[] data = new PositionInfo[]{
            new PositionInfo(43.48932915, 1.66561772, 1000, -20),
            new PositionInfo(43.48849093, 1.6648176, 2000, -10),
            new PositionInfo(43.48818612, 1.66615113, 3000, -50)
    };


    double[] target = new double[data.length];
    Arrays.fill(target, 0.0);

    double[] start = new double[2];

    for (PositionInfo row : data) {
        start[0] += row.latitude;
        start[1] += row.longitude;
    }
    start[0] /= data.length;
    start[1] /= data.length;

    MultivariateJacobianFunction distancesModel = new MultivariateJacobianFunction() {
        @Override
        public Pair<RealVector, RealMatrix> value(final RealVector point) {
            double tgtLat = point.getEntry(0);
            double tgtLong = point.getEntry(1);

            RealVector value = new ArrayRealVector(data.length);
            RealMatrix jacobian = new Array2DRowRealMatrix(data.length, 2);
            for (int i = 0; i < data.length; i++) {
                DerivativeStructure distance = getWeightedEuclideanDistance(tgtLat, tgtLong, data[i]);
                //DerivativeStructure distance = getWeightedHaversineDistance(tgtLat, tgtLong, data[i]);
                value.setEntry(i, distance.getValue());
                jacobian.setEntry(i, 0, distance.getPartialDerivative(1, 0));
                jacobian.setEntry(i, 1, distance.getPartialDerivative(0, 1));
            }

            return new Pair<RealVector, RealMatrix>(value, jacobian);
        }
    };


    LeastSquaresProblem problem = new LeastSquaresBuilder()
            .start(start)
            .model(distancesModel)
            .target(target)
            .lazyEvaluation(false)
            .maxEvaluations(1000)
            .maxIterations(1000)
            .build();

    LeastSquaresOptimizer optimizer = new LevenbergMarquardtOptimizer().
            withCostRelativeTolerance(1.0e-12).
            withParameterRelativeTolerance(1.0e-12);

    LeastSquaresOptimizer.Optimum optimum = optimizer.optimize(problem);
    RealVector point = optimum.getPoint();
    System.out.println("Start = " + Arrays.toString(start));
    System.out.println("Solve = " + point);
}

附注重量的逻辑对我来说似乎很可疑。在您引用的问题中，OP 对半径进行了一些估计，然后它是一个明显的权重。使用以对数测量的信号强度的倒平方 dBm我觉得很奇怪。