Matrix Square Root Example (omode01)
/***************************************************************/
/* */
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: omode01 */
/* TITLE: Matrix Square Root Example (omode01) */
/* PRODUCT: OR */
/* SYSTEM: ALL */
/* KEYS: OR */
/* PROCS: OPTMODEL */
/* DATA: */
/* */
/* SUPPORT: UPDATE: */
/* REF: */
/* MISC: Example 1 from the OPTMODEL chapter of */
/* Mathematical Programming. */
/* */
/***************************************************************/
proc optmodel;
number n = 5; /* size of matrix */
/* random original array */
number A{1..n, 1..n} = 10 - 20*rand('UNIFORM');
/* compute upper triangle of the
* symmetric matrix A*transpose(A) */
/* should be positive def unless A is singular */
number P{i in 1..n, j in i..n};
for {i in 1..n, j in i..n}
P[i,j] = sum{k in 1..n} A[i,k]*A[j,k];
/* coefficients of square root array
* (upper triangle of symmetric matrix) */
var q{i in 1..n, i..n};
/* The default initial value q[i,j]=0 is
* a local minimum of the objective,
* so you must move it away from that point. */
q[1,1] = 1;
/* minimize difference of square of q from P */
min r = sum{i in 1..n, j in i..n}
( sum{k in 1..i} q[k,i]*q[k,j]
+ sum{k in i+1..j} q[i,k]*q[k,j]
+ sum{k in j+1..n} q[i,k]*q[j,k]
- P[i,j] )**2;
solve;
print q;
quit;