r"""Compute the closest point (orthogonal projection) on the generalized `(n-1)`-simplex of a vector :math:`\mathbf{v}` wrt. to the Euclidean distance, thus solving: .. math:: \mathcal{P}(w) \in \mathop{\arg \min}_\gamma \| \gamma - \mathbf{v} \|_2 s.t. \ \gamma^T \math
(v, z=1)
| 84 | |
| 85 | |
| 86 | def proj_simplex(v, z=1): |
| 87 | r"""Compute the closest point (orthogonal projection) on the |
| 88 | generalized `(n-1)`-simplex of a vector :math:`\mathbf{v}` wrt. to the Euclidean |
| 89 | distance, thus solving: |
| 90 | |
| 91 | .. math:: |
| 92 | \mathcal{P}(w) \in \mathop{\arg \min}_\gamma \| \gamma - \mathbf{v} \|_2 |
| 93 | |
| 94 | s.t. \ \gamma^T \mathbf{1} = z |
| 95 | |
| 96 | \gamma \geq 0 |
| 97 | |
| 98 | If :math:`\mathbf{v}` is a 2d array, compute all the projections wrt. axis 0 |
| 99 | |
| 100 | .. note:: This function is backend-compatible and will work on arrays |
| 101 | from all compatible backends. |
| 102 | |
| 103 | Parameters |
| 104 | ---------- |
| 105 | v : {array-like}, shape (n, d) |
| 106 | z : int, optional |
| 107 | 'size' of the simplex (each vectors sum to z, 1 by default) |
| 108 | |
| 109 | Returns |
| 110 | ------- |
| 111 | h : ndarray, shape (`n`, `d`) |
| 112 | Array of projections on the simplex |
| 113 | """ |
| 114 | nx = get_backend(v) |
| 115 | n = v.shape[0] |
| 116 | if v.ndim == 1: |
| 117 | d1 = 1 |
| 118 | v = v[:, None] |
| 119 | else: |
| 120 | d1 = 0 |
| 121 | d = v.shape[1] |
| 122 | |
| 123 | # sort u in ascending order |
| 124 | u = nx.sort(v, axis=0) |
| 125 | # take the descending order |
| 126 | u = nx.flip(u, 0) |
| 127 | cssv = nx.cumsum(u, axis=0) - z |
| 128 | ind = nx.arange(n, type_as=v)[:, None] + 1 |
| 129 | cond = u - cssv / ind > 0 |
| 130 | rho = nx.sum(cond, 0) |
| 131 | theta = cssv[rho - 1, nx.arange(d)] / rho |
| 132 | w = nx.maximum(v - theta[None, :], nx.zeros(v.shape, type_as=v)) |
| 133 | if d1: |
| 134 | return w[:, 0] |
| 135 | else: |
| 136 | return w |
| 137 | |
| 138 | |
| 139 | def projection_sparse_simplex(V, max_nz, z=1, axis=None, nx=None): |