## Matrix Square Root (omod1)

```/***************************************************************/
/*                                                             */
/*          S A S   S A M P L E   L I B R A R Y                */
/*                                                             */
/*    NAME: omod1                                              */
/*   TITLE: Matrix Square Root (omod1)                         */
/* 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*ranuni(-1);
/* 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;

```