conjugate_gradient.sql_in

Go to the documentation of this file.

/* ----------------------------------------------------------------------- *//** 
 *
 * @file conjugate_gradient.sql_in
 *
 * @brief SQL function computing Conjugate Gradient
 * @date   March 2011
 *
 *
 *//* ----------------------------------------------------------------------- */

/**
@addtogroup grp_cg

@about
This function uses the iterative conjugate gradient method [1] to find a solution to the function: \f[ \boldsymbol Ax = \boldsymbol b \f]
where \f$ \boldsymbol A \f$ is a symmetric, positive definite matrix and \f$x\f$ and \f$ \boldsymbol b \f$ are vectors. 

@input
Matrix \f$ \boldsymbol A \f$ is assumed to be stored in a table where each row consists of at least two columns: array containing values of a given row, row number:
<pre>{TABLE|VIEW} <em>matrix_A</em> (
    <em>row_number</em> FLOAT,
    <em>row_values</em> FLOAT[],
)</pre>
The number of elements in each row should be the same.

\f$ \boldsymbol b \f$ is passed as a FLOAT[] to the function.

@usage
Conjugate gradient can be called as follows:
<pre>SELECT \ref conjugate_gradient('<em>table_name</em>', 
    '<em>name_of_row_values_col</em>', '<em>name_of_row_number_col</em>', '<em>aray_of_b_values</em>', 
    '<em>desired_precision</em>');</pre>
Function returns x as an array.  
    
@examp
-# Construct matrix A according to structure:
\code
sql> SELECT * FROM data;
 row_num | row_val 
---------+---------
       1 | {2,1}
       2 | {1,4}
(2 rows)
\endcode
-# Call conjugate gradient function:
\code
sql> SELECT conjugate_gradient('data','row_val','row_num','{2,1}',1E-6,1);
INFO:  COMPUTE RESIDUAL ERROR 14.5655661859659
INFO:  ERROR 0.144934004246004
INFO:  ERROR 3.12963615962926e-31
INFO:  TEST FINAL ERROR 2.90029642185163e-29
    conjugate_gradient     
---------------------------
 {1,-1.31838984174237e-15}
(1 row)
\endcode

@literature
[1] "Conjugate gradient method" Wikipedia - http://en.wikipedia.org/wiki/Conjugate_gradient_method

@sa File conjugate_gradient.sql_in documenting the SQL function.
*/

/**
 * @brief Compute conjugate gradient
 * 
 * @param matrix Name of the table containing argument matrix A
 * @param val_id Name of the column contains row values
 * @param row_id Name of the column contains row number
 * @param b Array containing values of b
 * @param precision_limit Precision threshold after which process will terminate
 * @param verbosity Verbose flag (0 = false, 1 = true)
 * @returns Array containing values of x
 *
 */
CREATE OR REPLACE FUNCTION MADLIB_SCHEMA.conjugate_gradient(Matrix TEXT, val_id TEXT, row_id TEXT, b FLOAT[], precision_limit FLOAT, verbosity INT)  RETURNS FLOAT[] AS $$
declare
    r FLOAT[];
    p FLOAT[];
    x FLOAT[];
    k INT;
    iter INT = 0;
    recidual_refresh INT := 30;
    alpha FLOAT;
    r_size FLOAT;
    r_new_size FLOAT;
    Ap FLOAT[];
    Ax FLOAT[];
    pAp_size FLOAT;
    beta FLOAT;
    exit_if_no_progess_in INT := 15;
begin   
    DROP TABLE IF EXISTS X_val;
    CREATE TEMP TABLE X_val(
        value FLOAT[]
    );
    
    DROP TABLE IF EXISTS P_val;
    CREATE TEMP TABLE P_val(
        value FLOAT[]
    ); 
    
    SELECT INTO k array_upper(b,1);
    --INSERT INTO X_val SELECT ARRAY(SELECT random() FROM generate_series(1, k));
    EXECUTE 'INSERT INTO X_val SELECT ARRAY(SELECT random() FROM generate_series(1,' ||k|| '))';
    LOOP
        IF(iter%recidual_refresh = 0)THEN 
            EXECUTE 'SELECT ARRAY(SELECT MADLIB_SCHEMA.array_dot('||val_id||', j.x) FROM (SELECT value AS x FROM X_val) AS j, '|| Matrix ||' ORDER BY '||row_id||' LIMIT '|| k ||')' INTO Ax;
            SELECT INTO r MADLIB_SCHEMA.array_sub(b, Ax);
            SELECT INTO r_size MADLIB_SCHEMA.array_dot(r, r);
            IF(verbosity > 0) THEN
                RAISE INFO 'COMPUTE RESIDUAL ERROR %', r_size;
            END IF;
            SELECT INTO p r; 
        END IF;
        iter = iter + 1;
        TRUNCATE TABLE P_val;
        --INSERT INTO P_val VALUES(p);
        EXECUTE 'INSERT INTO P_val(value) VALUES(array['|| array_to_string(p,',') ||'])';
        EXECUTE 'SELECT ARRAY(SELECT MADLIB_SCHEMA.array_dot('||val_id||', j.p) FROM (SELECT value AS p FROM P_val) AS j,'|| Matrix ||' ORDER BY '||row_id||' LIMIT '|| k ||')' INTO Ap;
        SELECT INTO pAp_size MADLIB_SCHEMA.array_dot(p, Ap);
        alpha = r_size/pAp_size;
        
        --SELECT INTO x MADLIB_SCHEMA.array_add(value, MADLIB_SCHEMA.array_scalar_mult(p,alpha)) FROM X_val;
        EXECUTE 'SELECT MADLIB_SCHEMA.array_add(value, MADLIB_SCHEMA.array_scalar_mult(array['|| array_to_string(p,',') ||']::float[],'||alpha||'::float)) FROM X_val' INTO x;
        TRUNCATE TABLE X_val;
        --INSERT INTO X_val VALUES(x);
        EXECUTE 'INSERT INTO X_val VALUES(array['|| array_to_string(x,',') ||'])';
        
        SELECT INTO r MADLIB_SCHEMA.array_add(r,MADLIB_SCHEMA.array_scalar_mult(Ap, -alpha));
        SELECT INTO r_new_size MADLIB_SCHEMA.array_dot(r,r);
        
        IF(verbosity > 0) THEN
            RAISE INFO 'ERROR %',r_new_size; 
        END IF;
        IF (r_new_size < precision_limit) THEN
            EXECUTE 'SELECT ARRAY(SELECT MADLIB_SCHEMA.array_dot('||val_id||', j.x) FROM (SELECT value AS x FROM X_val) AS j, '|| Matrix ||' ORDER BY '||row_id||' LIMIT '|| k ||')' INTO Ax;
            SELECT INTO r MADLIB_SCHEMA.array_sub(b, Ax);
            SELECT INTO r_new_size MADLIB_SCHEMA.array_dot(r, r);
            IF(verbosity > 0) THEN
                RAISE INFO 'TEST FINAL ERROR %', r_new_size;
            END IF;
            IF (r_new_size < precision_limit) THEN
                EXIT;
            END IF;
        END IF;
        SELECT INTO p MADLIB_SCHEMA.array_add(r, MADLIB_SCHEMA.array_scalar_mult(p, r_new_size/r_size));
        IF(r_size < r_new_size) THEN
            exit_if_no_progess_in = exit_if_no_progess_in-1;
            RAISE INFO 'No progress! count = %',exit_if_no_progess_in;
            
            IF(exit_if_no_progess_in <= 0) THEN
                RAISE EXCEPTION 'Algorithm failed to converge. Check is input is positive definate.';
            END IF;
        ELSE
            exit_if_no_progess_in = 15;
        END IF;
        r_size = r_new_size;
    END LOOP; 
    IF(verbosity > 1) THEN
        RETURN ARRAY[r_new_size];
    END IF;
    --SELECT INTO x value FROM X_val;
    EXECUTE 'SELECT value FROM X_val' INTO x;
    RETURN x;
end
$$ LANGUAGE plpgsql;

CREATE OR REPLACE FUNCTION MADLIB_SCHEMA.conjugate_gradient(Matrix TEXT, val_id TEXT, row_id TEXT, b FLOAT[], precision_limit FLOAT)  RETURNS FLOAT[] AS $$
declare
begin
    RETURN MADLIB_SCHEMA.conjugate_gradient(Matrix, val_id, row_id, b, precision_limit,0);
end
$$ LANGUAGE plpgsql;