Commit d92cb1e6 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

spectral: add Bethe Hessian implementation

parent bca4052a
......@@ -32,7 +32,7 @@ using namespace boost;
using namespace graph_tool;
void laplacian(GraphInterface& g, boost::any index, boost::any weight,
string sdeg,
string sdeg, double r,
python::object odata, python::object oi,
python::object oj)
{
......@@ -67,13 +67,14 @@ void laplacian(GraphInterface& g, boost::any index, boost::any weight,
return get_laplacian()
(std::forward<decltype(graph)>(graph),
std::forward<decltype(a2)>(a2),
std::forward<decltype(a3)>(a3), deg, data, i, j);
std::forward<decltype(a3)>(a3), deg, r, data, i, j);
},
vertex_scalar_properties(), weight_props_t())(index, weight);
}
void laplacian_matvec(GraphInterface& g, boost::any index, boost::any weight,
boost::any deg, python::object ov, python::object oret)
boost::any deg, double r, python::object ov,
python::object oret)
{
if (!belongs<vertex_scalar_properties>()(index))
throw ValueException("index vertex property must have a scalar value type");
......@@ -98,13 +99,14 @@ void laplacian_matvec(GraphInterface& g, boost::any index, boost::any weight,
(g,
[&](auto&& graph, auto&& vi, auto&& w)
{
return lap_matvec(graph, vi, w, d, v, ret);
return lap_matvec(graph, vi, w, d, r, v, ret);
},
vertex_scalar_properties(), weight_props_t())(index, weight);
}
void laplacian_matmat(GraphInterface& g, boost::any index, boost::any weight,
boost::any deg, python::object ov, python::object oret)
boost::any deg, double r, python::object ov,
python::object oret)
{
if (!belongs<vertex_scalar_properties>()(index))
throw ValueException("index vertex property must have a scalar value type");
......@@ -129,7 +131,7 @@ void laplacian_matmat(GraphInterface& g, boost::any index, boost::any weight,
(g,
[&](auto&& graph, auto&& vi, auto&& w)
{
return lap_matmat(graph, vi, w, d, v, ret);
return lap_matmat(graph, vi, w, d, r, v, ret);
},
vertex_scalar_properties(), weight_props_t())(index, weight);
}
......@@ -73,6 +73,7 @@ struct get_laplacian
{
template <class Graph, class Index, class Weight>
void operator()(const Graph& g, Index index, Weight weight, deg_t deg,
double r,
multi_array_ref<double,1>& data,
multi_array_ref<int32_t,1>& i,
multi_array_ref<int32_t,1>& j) const
......@@ -83,20 +84,21 @@ struct get_laplacian
if (source(e, g) == target(e, g))
continue;
data[pos] = -get(weight, e);
data[pos] = -get(weight, e) * r;
i[pos] = get(index, target(e, g));
j[pos] = get(index, source(e, g));
++pos;
if (!graph_tool::is_directed(g))
{
data[pos] = -get(weight, e);
data[pos] = -get(weight, e) * r;
i[pos] = get(index, source(e, g));
j[pos] = get(index, target(e, g));
++pos;
}
}
double rr = r*r - 1;
for (auto v : vertices_range(g))
{
double k = 0;
......@@ -111,7 +113,7 @@ struct get_laplacian
case TOTAL_DEG:
k = sum_degree(g, v, weight, all_edges_iteratorS<Graph>());
}
data[pos] = k;
data[pos] = k + rr;
i[pos] = j[pos] = get(index, v);
++pos;
}
......@@ -120,8 +122,9 @@ struct get_laplacian
};
template <class Graph, class Vindex, class Weight, class Deg, class V>
void lap_matvec(Graph& g, Vindex index, Weight w, Deg d, V& x, V& ret)
void lap_matvec(Graph& g, Vindex index, Weight w, Deg d, double r, V& x, V& ret)
{
double rr = r*r - 1;
parallel_vertex_loop
(g,
[&](auto v)
......@@ -135,7 +138,7 @@ void lap_matvec(Graph& g, Vindex index, Weight w, Deg d, V& x, V& ret)
if (u == v)
continue;
auto w_e = get(w, e);
y += w_e * x[get(index, u)];
y += r * w_e * x[get(index, u)];
}
}
else
......@@ -144,16 +147,17 @@ void lap_matvec(Graph& g, Vindex index, Weight w, Deg d, V& x, V& ret)
{
if (u == v)
continue;
y += x[get(index, u)];
y += r * x[get(index, u)];
}
}
ret[get(index, v)] = d[v] * x[get(index, v)] - y;
ret[get(index, v)] = (d[v] + rr) * x[get(index, v)] - y;
});
}
template <class Graph, class Vindex, class Weight, class Deg, class Mat>
void lap_matmat(Graph& g, Vindex index, Weight w, Deg d, Mat& x, Mat& ret)
void lap_matmat(Graph& g, Vindex index, Weight w, Deg d, double r, Mat& x, Mat& ret)
{
double rr = r*r - 1;
size_t M = x.shape()[1];
parallel_vertex_loop
(g,
......@@ -171,7 +175,7 @@ void lap_matmat(Graph& g, Vindex index, Weight w, Deg d, Mat& x, Mat& ret)
auto w_e = get(w, e);
auto ui = get(index, u);
for (size_t i = 0; i < M; ++i)
y[i] += w_e * x[ui][i];
y[i] += r * w_e * x[ui][i];
}
}
else
......@@ -182,11 +186,11 @@ void lap_matmat(Graph& g, Vindex index, Weight w, Deg d, Mat& x, Mat& ret)
continue;
auto ui = get(index, u);
for (size_t i = 0; i < M; ++i)
y[i] += x[ui][i];
y[i] += r * x[ui][i];
}
}
for (size_t i = 0; i < M; ++i)
ret[vi][i] = d[v] * x[vi][i] - y[i];
ret[vi][i] = (d[v] + rr) * x[vi][i] - y[i];
});
}
......
......@@ -33,15 +33,17 @@ void adjacency_matmat(GraphInterface& g, boost::any index, boost::any weight,
python::object ov, python::object oret);
void laplacian(GraphInterface& g, boost::any index, boost::any weight,
string sdeg,
string sdeg, double r,
python::object odata, python::object oi,
python::object oj);
void laplacian_matvec(GraphInterface& g, boost::any index, boost::any weight,
boost::any deg, python::object ov, python::object oret);
boost::any deg, double r, python::object ov,
python::object oret);
void laplacian_matmat(GraphInterface& g, boost::any index, boost::any weight,
boost::any deg, python::object ov, python::object oret);
boost::any deg, double r, python::object ov,
python::object oret);
void norm_laplacian(GraphInterface& g, boost::any index, boost::any weight,
string sdeg,
......
......@@ -252,9 +252,9 @@ class AdjacencyOperator(scipy.sparse.linalg.LinearOperator):
return self
@_limit_args({"deg": ["total", "in", "out"]})
def laplacian(g, deg="out", norm=False, weight=None, vindex=None, operator=False,
def laplacian(g, deg="out", norm=False, weight=None, r=1, vindex=None, operator=False,
csr=True):
r"""Return the Laplacian matrix of the graph.
r"""Return the Laplacian (or Bethe Hessian if :math:`r > 1`) matrix of the graph.
Parameters
----------
......@@ -266,6 +266,10 @@ def laplacian(g, deg="out", norm=False, weight=None, vindex=None, operator=False
Whether to compute the normalized Laplacian.
weight : :class:`~graph_tool.EdgePropertyMap` (optional, default: ``None``)
Edge property map with the edge weights.
r : ``double`` (optional, default: ``1.``)
Regularization parameter. If :math:`r > 1`, and ``norm`` is ``False``,
then this corresponds to the Bethe Hessian. (This parameter has an
effect only for ``norm == False``.)
vindex : :class:`~graph_tool.VertexPropertyMap` (optional, default: ``None``)
Vertex property map specifying the row/column indexes. If not provided, the
internal vertex index is used.
......@@ -296,7 +300,21 @@ def laplacian(g, deg="out", norm=False, weight=None, vindex=None, operator=False
\end{cases}
Where :math:`\Gamma(v_i)=\sum_j A_{ij}w_{ij}` is sum of the weights of the
edges incident on vertex :math:`v_i`. The normalized version is
edges incident on vertex :math:`v_i`.
In case of :math:`r > 1`, the matrix returned is the Bethe Hessian
[bethe-hessian]_:
.. math::
\ell_{ij} =
\begin{cases}
\Gamma(v_i) + (r^2 - 1) & \text{if } i = j \\
-r w_{ij} & \text{if } i \neq j \text{ and } (j, i) \in E \\
0 & \text{otherwise}.
\end{cases}
The normalized version is
.. math::
......@@ -386,10 +404,15 @@ def laplacian(g, deg="out", norm=False, weight=None, vindex=None, operator=False
References
----------
.. [wikipedia-laplacian] http://en.wikipedia.org/wiki/Laplacian_matrix
.. [bethe-hessian] Saade, Alaa, Florent Krzakala, and Lenka
Zdeborová. "Spectral clustering of graphs with the bethe hessian." Advances
in Neural Information Processing Systems 27 (2014): 406-414, :arxiv:`1406.1880`,
https://proceedings.neurips.cc/paper/2014/hash/63923f49e5241343aa7acb6a06a751e7-Abstract.html
"""
if operator:
return LaplacianOperator(g, deg=deg, norm=norm, weight=weight,
return LaplacianOperator(g, deg=deg, norm=norm, weight=weight, r=r,
vindex=vindex)
if vindex is None:
......@@ -415,7 +438,7 @@ def laplacian(g, deg="out", norm=False, weight=None, vindex=None, operator=False
_prop("e", g, weight), deg, data, i, j)
else:
libgraph_tool_spectral.laplacian(g._Graph__graph, _prop("v", g, vindex),
_prop("e", g, weight), deg, data, i, j)
_prop("e", g, weight), deg, r, data, i, j)
if E > 0:
V = max(g.num_vertices(), max(i.max() + 1, j.max() + 1))
else:
......@@ -429,12 +452,13 @@ def laplacian(g, deg="out", norm=False, weight=None, vindex=None, operator=False
class LaplacianOperator(scipy.sparse.linalg.LinearOperator):
@_limit_args({"deg": ["total", "in", "out"]})
def __init__(self, g, weight=None, deg="out", norm=False, vindex=None):
def __init__(self, g, weight=None, deg="out", r=1, norm=False, vindex=None):
r"""A :class:`scipy.sparse.linalg.LinearOperator` representing the laplacian
matrix of a graph. See :func:`laplacian` for details."""
self.g = g
self.weight = weight
self.r = r
if vindex is None:
if g.get_vertex_filter()[0] is not None:
......@@ -469,12 +493,16 @@ class LaplacianOperator(scipy.sparse.linalg.LinearOperator):
x = numpy.asarray(x, dtype="float")
if self.norm:
matvec = libgraph_tool_spectral.norm_laplacian_matvec
matvec(self.g._Graph__graph,
_prop("v", self.g, self.vindex),
_prop("e", self.g, self.weight),
_prop("v", self.g, self.d), x, y)
else:
matvec = libgraph_tool_spectral.laplacian_matvec
matvec(self.g._Graph__graph,
_prop("v", self.g, self.vindex),
_prop("e", self.g, self.weight),
_prop("v", self.g, self.d), x, y)
matvec(self.g._Graph__graph,
_prop("v", self.g, self.vindex),
_prop("e", self.g, self.weight),
_prop("v", self.g, self.d), self.r, x, y)
return y
def _matmat(self, x):
......@@ -482,12 +510,16 @@ class LaplacianOperator(scipy.sparse.linalg.LinearOperator):
x = numpy.asarray(x, dtype="float")
if self.norm:
matmat = libgraph_tool_spectral.norm_laplacian_matmat
matmat(self.g._Graph__graph,
_prop("v", self.g, self.vindex),
_prop("e", self.g, self.weight),
_prop("v", self.g, self.d), x, y)
else:
matmat = libgraph_tool_spectral.laplacian_matmat
matmat(self.g._Graph__graph,
_prop("v", self.g, self.vindex),
_prop("e", self.g, self.weight),
_prop("v", self.g, self.d), x, y)
matmat(self.g._Graph__graph,
_prop("v", self.g, self.vindex),
_prop("e", self.g, self.weight),
_prop("v", self.g, self.d), self.r, x, y)
return y
def _adjoint(self):
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment