Fit the D1 of a MAP by MINIMIZING a distance to a reference MAP. More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <vector>
#include "line/api/mam/map_dist.h"
#include "line/api/mam/map_moment.h"
#include "line/util/lu.h"
#include "line/num/number.h"
#include "line/util/auglag.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
#include "line/util/sylvester.h"

Include dependency graph for map_optim_dist.h:

Classes
struct	line::mam::MapOptimDist< T >
	What an optimizing distance fit returns. More...

Namespaces
namespace	line
namespace	line::mam

Functions
template<class T>
MapOptimDist< T >	line::mam::map_optim_dist (const Map< T > &a, const std::vector< T > &alA, const Matrix< T > &B0, const std::vector< T > &alB, unsigned L)
	Fit B1 minimizing the lag-L joint-density distance to a, with B0 fixed.
template<class T>
MapOptimDist< T >	line::mam::map_optim_dist_acf (const Map< T > &a, const std::vector< T > &alA, const Matrix< T > &B0, const std::vector< T > &alB)
	Fit B1 minimizing the AUTOCORRELATION distance to a, with B0 fixed.

Detailed Description

Fit the D1 of a MAP by MINIMIZING a distance to a reference MAP.

Port of matlab/lib/kpctoolbox/map/map_optim_dist.m and map_optim_dist_acf.m (twins in python/line_solver/api/mapdist/continuous.py). The distances themselves – map_dist, map_dist_acf – are already in map_dist.h; what is here is the OPTIMIZATION over D1 with D0 held fixed.

Reference: G. Horvath, "Measuring the distance between MAPs and some applications", ASMTA 2015, LNCS 9081, pp. 95-109.

WHAT IS BEING FITTED, AND WHY D0 IS FIXED. Given a reference MAP A and a chosen D0 for the approximation B, the free parameter is B's D1. The constraints are exactly the two that make (B0, B1) a MAP with the DECLARED embedded distribution alB:

alB (-B0)^-1 B1 = alB, i.e. alB is stationary at arrivals. This is the kron(I, alB inv(-B0)) block;
each row of B0 + B1 sums to zero, i.e. B1's row sums are -B0's. This is the kron(ones(1,N), I) block.

Every entry is bounded below by 1e-6 rather than by 0, which is the reference's own floor: an exactly zero entry makes the fitted MAP reducible, and the distance is then defined on a different state space than the one the constraints were written for.

THE LAG-1 CASE IS A CONVEX QP AND IS SOLVED AS ONE. At L = 1 the distance is a quadratic form in vec(B1),

d(vB1) = vB1' H vB1 + vA1' H_AA vA1 - 2 vA1' H_AB vB1,

with H = mkron(X_BB, Z_BB) positive semidefinite, so the reference calls quadprog and gets a global optimum. Every other case is a general nonlinear program over the same affine set and the reference calls fmincon, which is local. Both are served here by auglag, and the DISTINCTION IS PRESERVED IN WHAT THE RESULT PROMISES: global is set only on the L = 1 path.

The six Lyapunov solves are the reference's; lyap_solve is the tree's MATLAB-compatible lyap(A,B,C), solving A X + X B + C = 0.

COLUMN-MAJOR THROUGHOUT. vec here is MATLAB's, stacking COLUMNS, because the Kronecker identity the quadratic form rests on (vec(A X B) = kron(B', A) vec(X)) holds in that convention and in no other. Getting this wrong transposes the fitted D1 without changing its row sums, so the constraints still pass and only the distance is wrong.

ARITHMETIC: transcendental, inherited from the distances.

Definition in file map_optim_dist.h.

Classes

Namespaces

Functions

Detailed Description