Skip to content

Commit

Permalink
brain fart
Browse files Browse the repository at this point in the history
  • Loading branch information
araujoms committed Dec 28, 2024
1 parent 584f88e commit d43d1ce
Show file tree
Hide file tree
Showing 2 changed files with 17 additions and 17 deletions.
22 changes: 11 additions & 11 deletions docs/src/api/cones.rst
Original file line number Diff line number Diff line change
Expand Up @@ -88,9 +88,9 @@ and for short we use :math:`S \succeq 0` to denote membership. SCS
vectorizes these cones in a special way which we detail here.

SCS assumes that the input data corresponding to positive semidefinite cones have been vectorized by scaling the off-diagonal entries by
:math:`\sqrt{2}` and stacking the real parameters of the upper triangular row-wise. For a :math:`k \times k`
real matrix variable (or data matrix) this operation would create a vector of length
:math:`k(k+1)/2`, whereas for a :math:`k \times k` complex matrix its creates a vector of length :math:`k^2`. Scaling by :math:`\sqrt{2}` is required to preserve the inner-product.
:math:`\sqrt{2}` and stacking the real parameters of the lower triangle column-wise. For a :math:`k \times k`
real matrix variable (or data matrix) this operation creates a vector of length
:math:`k(k+1)/2`, whereas for a :math:`k \times k` complex matrix it creates a vector of length :math:`k^2`. Scaling by :math:`\sqrt{2}` is required to preserve the inner-product.

**This must be done for the rows of both** :math:`A` **and** :math:`b` **that correspond to semidefinite cones and must be done independently for each semidefinite cone.**

Expand All @@ -105,20 +105,20 @@ where the :math:`\text{vec}` operation takes the :math:`k \times k` matrix
S_{k1} & S_{k2} & \ldots & S_{kk} \\
\end{bmatrix}
and produces a vector consisting of the upper triangular real parameters elements scaled and re-arranged. In the real case, the result is
and produces a vector consisting of the lower triangular real parameters scaled and re-arranged. In the real case, the result is

.. math::
\text{vec}(S) = (S_{11}, \sqrt{2} S_{12}, \ldots, \sqrt{2} S_{1k}, S_{22}, \sqrt{2}S_{23}, \dots, S_{k-1,k-1}, \sqrt{2}S_{k-1,k}, S_{kk}) \in \mathbf{R}^{k(k+1)/2},
\text{vec}(S) = (S_{11}, \sqrt{2} S_{21}, \ldots, \sqrt{2} S_{k1}, S_{22}, \sqrt{2}S_{32}, \dots, S_{k-1,k-1}, \sqrt{2}S_{k,k-1}, S_{kk}) \in \mathbf{R}^{k(k+1)/2},
whereas in the complex case the result is

.. math::
\text{vec}(S) = (S_{11}, \sqrt{2} \Re(S_{12}), \sqrt{2} \Im(S_{12}), \ldots, \sqrt{2} \Re(S_{1k}), \sqrt{2} \Im(S_{1k}), S_{22}, \\\sqrt{2}\Re(S_{23}), \sqrt{2}\Im(S_{23}), \dots, S_{k-1,k-1}, \sqrt{2}\Re(S_{k-1,k}), \sqrt{2}\Im(S_{k-1,k}), S_{kk}) \in \mathbf{R}^{k^2}.
\text{vec}(S) = (S_{11}, \sqrt{2} \Re(S_{21}), \sqrt{2} \Im(S_{21}), \ldots, \sqrt{2} \Re(S_{k1}), \sqrt{2} \Im(S_{k1}), S_{22}, \\\sqrt{2}\Re(S_{32}), \sqrt{2}\Im(S_{32}), \dots, S_{k-1,k-1}, \sqrt{2}\Re(S_{k,k-1}), \sqrt{2}\Im(S_{k,k-1}), S_{kk}) \in \mathbf{R}^{k^2}.
To recover the matrix solution this operation must be inverted on the components
of the vectors returned by SCS corresponding to each semidefinite cone. That is, the
off-diagonal entries must be scaled by :math:`1/\sqrt{2}` and the lower triangular
entries are filled in by copying the values of upper triangular entries.
off-diagonal entries must be scaled by :math:`1/\sqrt{2}` and the upper triangular
entries are filled in by copying the values of lower triangular entries.
Explicitly, in the real case the inverse operation takes vector :math:`s \in
\mathbf{R}^{k(k+1)/2}` and produces the matrix

Expand All @@ -136,10 +136,10 @@ whereas in the complex case the inverse operation takes vector :math:`s \in

.. math::
\text{mat}(s) = \begin{bmatrix}
s_{1} & (s_{2} + i s_3) / \sqrt{2} & \ldots & (s_{2k-2}+is_{2k-1}) / \sqrt{2} \\
(s_{2} - i s_3) / \sqrt{2} / \sqrt{2} & s_{2k} & \ldots & (s_{4k-5}+is_{4k-4}) / \sqrt{2} \\
s_{1} & (s_{2} - i s_3) / \sqrt{2} & \ldots & (s_{2k-2} - is_{2k-1}) / \sqrt{2} \\
(s_{2} + i s_3) / \sqrt{2} & s_{2k} & \ldots & (s_{4k-5} - is_{4k-4}) / \sqrt{2} \\
\vdots & \vdots & \ddots & \vdots \\
(s_{2k-2}-is_{2k-1}) / \sqrt{2} & (s_{4k-5}-is_{4k-4}) / \sqrt{2} & \ldots & s_{k^2} \\
(s_{2k-2} + is_{2k-1}) / \sqrt{2} & (s_{4k-5} + is_{4k-4}) / \sqrt{2} & \ldots & s_{k^2} \\
\end{bmatrix}
\in \mathbf{C}^{k \times k},
Expand Down
12 changes: 6 additions & 6 deletions src/cones.c
Original file line number Diff line number Diff line change
Expand Up @@ -487,7 +487,7 @@ static scs_int proj_semi_definite_cone(scs_float *X, const scs_int n,

#ifdef USE_LAPACK

/* copy upper triangular matrix into full matrix */
/* copy lower triangular matrix into full matrix */
for (i = 0; i < n; ++i) {
memcpy(&(Xs[i * (n + 1)]), &(X[i * n - ((i - 1) * i) / 2]),
(n - i) * sizeof(scs_float));
Expand Down Expand Up @@ -539,7 +539,7 @@ static scs_int proj_semi_definite_cone(scs_float *X, const scs_int n,
/* undo rescaling: scale diags by 1/sqrt(2) */
BLAS(scal)(&nb, &sqrt2_inv, Xs, &nb_plus_one); /* not n_squared */

/* extract just upper triangular matrix */
/* extract just lower triangular matrix */
for (i = 0; i < n; ++i) {
memcpy(&(X[i * n - ((i - 1) * i) / 2]), &(Xs[i * (n + 1)]),
(n - i) * sizeof(scs_float));
Expand Down Expand Up @@ -648,7 +648,7 @@ static scs_int proj_complex_semi_definite_cone(scs_float *X, const scs_int n,

#ifdef USE_LAPACK

/* copy upper triangular matrix into full matrix */
/* copy lower triangular matrix into full matrix */
for (i = 0; i < n - 1; ++i) {
cXs[i * (n + 1)] = X[i * (2 * n - i)];
memcpy(&(cXs[i * (n + 1) + 1]), &(X[i * (2 * n - i) + 1]),
Expand Down Expand Up @@ -687,21 +687,21 @@ static scs_int proj_complex_semi_definite_cone(scs_float *X, const scs_int n,
return 0;
}

/*LAPACK is col-major, so the columns of cZ' are the eigenvectors */
/* cZ is matrix of all eigenvectors */
/* scale cZ by sqrt(eig) */
for (i = first_idx; i < n; ++i) {
csq_eig_pos = SQRTF(e[i]);
BLASC(scal)(&nb, &csq_eig_pos, &cZ[i * n], &one_int);
}

/* Xs = Z Z' = V E V' */
/* Xs = cZ cZ' = V E V' */
ncols_z = (blas_int)(n - first_idx);
BLASC(herk)("Lower", "NoTrans", &nb, &ncols_z, &one, &cZ[first_idx * n], &nb, &zero, cXs, &nb);

/* undo rescaling: scale diags by 1/sqrt(2) */
BLASC(scal)(&nb, &csqrt2_inv, cXs, &nb_plus_one); /* not n_squared */

/* extract just upper triangular matrix */
/* extract just lower triangular matrix */
for (i = 0; i < n - 1; ++i) {
X[i * (2 * n - i)] = cXs[i * (n + 1)];
memcpy(&(X[i * (2 * n - i) + 1]), &(cXs[i * (n + 1) + 1]),
Expand Down

0 comments on commit d43d1ce

Please sign in to comment.