Python scipy.optimize.leastsq для Java org.apache.commons.math3.fitting.leastsquares

Question

Я пытаюсь имитировать этот алгоритм, разработанный в Python, который рассчитывает геолокацию на основе местоположений видимых станций Wi-Fi, сам на основе этой идеи.

Этот алгоритм сначала использует функцию Numpy для вычисления среднего взвешенного значения наблюдаемых широт и долгот. Чтобы минимизировать влияние возможных ошибок в позициях Wi-Fi, он также использует метод " scipy.optimize.leastsq", чтобы рассчитать статистическим способом и, если возможно, более точную позицию.

Я хочу реализовать то же поведение на платформе Java Android.

Для всех остальных расчетов я успешно использую org.apache.commons.math3. Поэтому для решения задачи наименьших квадратов я логически пытаюсь опираться на https://commons.apache.org/proper/commons-math/userguide/leastsquares.html.

Моя проблема, если я хорошо понимаю, заключается в том, что Сципи управляет для меня сложностью определения функции Якоби, а мои плохие математические навыки не позволяют мне правильно определить модель LeastSquaresProblem. Я попробовал некоторые эксперименты, основанные на этом примере, которые кажутся закрытыми для того, что мне нужно, но результаты не очень хорошие, так как я не знаю, как обращаться с "якобианскими" частями.

Как кто-то делает для этого поста, может ли кто-то сделать то же самое для меня и попытаться объяснить это простым способом?

Подробнее о том, как работает 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 = 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, но я не уверен, что правильно его использую, так как не понимаю, как это работает.

К вашему сведению, я использую расчет расстояния Винсенти, который, например, можно заменить на Haversine или Euclidean.

Спасибо за любую помощь!

Answer 1

Код портировать нелегко, потому что SciPy предоставляет более общий интерфейс минимизации наименьших квадратов, а Apache Commons Math обеспечивает подгонку кривой. Тем не менее, многие проблемы оптимизации могут быть сформулированы как подгонка кривой. В коде 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 явно требует, чтобы вы предоставили оптимизатору средства для вычисления производных оптимизируемой функции. Код Python, похоже, использует какой-то дифференциатор различий, если Dfun не предоставляется Вы можете сделать что-то подобное в 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);
}

PS логика веса мне кажется подозрительной. В вопросе, на который вы ссылаетесь, ОП имеет некоторые оценки радиуса, а затем это очевидный вес. Использование обратного квадрата силы сигнала, которое измеряется в логарифмическом дБм, мне кажется странным.