0001: /*
0002: * This program is free software; you can redistribute it and/or modify
0003: * it under the terms of the GNU General Public License as published by
0004: * the Free Software Foundation; either version 2 of the License, or (at
0005: * your option) any later version.
0006: *
0007: * This program is distributed in the hope that it will be useful, but
0008: * WITHOUT ANY WARRANTY; without even the implied warranty of
0009: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
0010: * General Public License for more details.
0011: *
0012: * You should have received a copy of the GNU General Public License
0013: * along with this program; if not, write to the Free Software
0014: * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
0015:
0016: /*
0017: * PaceMatrix.java
0018: * Copyright (C) 2002 University of Waikato, Hamilton, New Zealand
0019: *
0020: */
0021:
0022: package weka.classifiers.functions.pace;
0023:
0024: import weka.core.matrix.DoubleVector;
0025: import weka.core.matrix.FlexibleDecimalFormat;
0026: import weka.core.matrix.IntVector;
0027: import weka.core.matrix.Matrix;
0028: import weka.core.matrix.Maths;
0029:
0030: import java.util.Random;
0031: import java.text.DecimalFormat;
0032:
0033: /**
0034: * Class for matrix manipulation used for pace regression. <p>
0035: *
0036: * REFERENCES <p>
0037: *
0038: * Wang, Y. (2000). "A new approach to fitting linear models in high
0039: * dimensional spaces." PhD Thesis. Department of Computer Science,
0040: * University of Waikato, New Zealand. <p>
0041: *
0042: * Wang, Y. and Witten, I. H. (2002). "Modeling for optimal probability
0043: * prediction." Proceedings of ICML'2002. Sydney. <p>
0044: *
0045: * @author Yong Wang (yongwang@cs.waikato.ac.nz)
0046: * @version $Revision: 1.5 $
0047: */
0048: public class PaceMatrix extends Matrix {
0049:
0050: /** for serialization */
0051: static final long serialVersionUID = 2699925616857843973L;
0052:
0053: /* ------------------------
0054: Constructors
0055: * ------------------------ */
0056:
0057: /** Construct an m-by-n PACE matrix of zeros.
0058: @param m Number of rows.
0059: @param n Number of colums.
0060: */
0061: public PaceMatrix(int m, int n) {
0062: super (m, n);
0063: }
0064:
0065: /** Construct an m-by-n constant PACE matrix.
0066: @param m Number of rows.
0067: @param n Number of colums.
0068: @param s Fill the matrix with this scalar value.
0069: */
0070: public PaceMatrix(int m, int n, double s) {
0071: super (m, n, s);
0072: }
0073:
0074: /** Construct a PACE matrix from a 2-D array.
0075: @param A Two-dimensional array of doubles.
0076: @throws IllegalArgumentException All rows must have the same length
0077: */
0078: public PaceMatrix(double[][] A) {
0079: super (A);
0080: }
0081:
0082: /** Construct a PACE matrix quickly without checking arguments.
0083: @param A Two-dimensional array of doubles.
0084: @param m Number of rows.
0085: @param n Number of colums.
0086: */
0087: public PaceMatrix(double[][] A, int m, int n) {
0088: super (A, m, n);
0089: }
0090:
0091: /** Construct a PaceMatrix from a one-dimensional packed array
0092: @param vals One-dimensional array of doubles, packed by columns (ala Fortran).
0093: @param m Number of rows.
0094: @throws IllegalArgumentException Array length must be a multiple of m.
0095: */
0096: public PaceMatrix(double vals[], int m) {
0097: super (vals, m);
0098: }
0099:
0100: /** Construct a PaceMatrix with a single column from a DoubleVector
0101: @param v DoubleVector
0102: */
0103: public PaceMatrix(DoubleVector v) {
0104: this (v.size(), 1);
0105: setMatrix(0, v.size() - 1, 0, v);
0106: }
0107:
0108: /** Construct a PaceMatrix from a Matrix
0109: @param X Matrix
0110: */
0111: public PaceMatrix(Matrix X) {
0112: super (X.getRowDimension(), X.getColumnDimension());
0113: A = X.getArray();
0114: }
0115:
0116: /* ------------------------
0117: Public Methods
0118: * ------------------------ */
0119:
0120: /** Set the row dimenion of the matrix
0121: * @param rowDimension the row dimension
0122: */
0123: public void setRowDimension(int rowDimension) {
0124: m = rowDimension;
0125: }
0126:
0127: /** Set the column dimenion of the matrix
0128: * @param columnDimension the column dimension
0129: */
0130: public void setColumnDimension(int columnDimension) {
0131: n = columnDimension;
0132: }
0133:
0134: /**
0135: * Clone the PaceMatrix object.
0136: *
0137: * @return the clone
0138: */
0139: public Object clone() {
0140: PaceMatrix X = new PaceMatrix(m, n);
0141: double[][] C = X.getArray();
0142: for (int i = 0; i < m; i++) {
0143: for (int j = 0; j < n; j++) {
0144: C[i][j] = A[i][j];
0145: }
0146: }
0147: return (Object) X;
0148: }
0149:
0150: /** Add a value to an element and reset the element
0151: * @param i the row number of the element
0152: * @param j the column number of the element
0153: * @param s the double value to be added with
0154: */
0155: public void setPlus(int i, int j, double s) {
0156: A[i][j] += s;
0157: }
0158:
0159: /** Multiply a value with an element and reset the element
0160: * @param i the row number of the element
0161: * @param j the column number of the element
0162: * @param s the double value to be multiplied with
0163: */
0164: public void setTimes(int i, int j, double s) {
0165: A[i][j] *= s;
0166: }
0167:
0168: /** Set the submatrix A[i0:i1][j0:j1] with a same value
0169: * @param i0 the index of the first element of the column
0170: * @param i1 the index of the last element of the column
0171: * @param j0 the index of the first column
0172: * @param j1 the index of the last column
0173: * @param s the value to be set to
0174: */
0175: public void setMatrix(int i0, int i1, int j0, int j1, double s) {
0176: try {
0177: for (int i = i0; i <= i1; i++) {
0178: for (int j = j0; j <= j1; j++) {
0179: A[i][j] = s;
0180: }
0181: }
0182: } catch (ArrayIndexOutOfBoundsException e) {
0183: throw new ArrayIndexOutOfBoundsException(
0184: "Index out of bounds");
0185: }
0186: }
0187:
0188: /** Set the submatrix A[i0:i1][j] with the values stored in a
0189: * DoubleVector
0190: * @param i0 the index of the first element of the column
0191: * @param i1 the index of the last element of the column
0192: * @param j the index of the column
0193: * @param v the vector that stores the values*/
0194: public void setMatrix(int i0, int i1, int j, DoubleVector v) {
0195: for (int i = i0; i <= i1; i++) {
0196: A[i][j] = v.get(i - i0);
0197: }
0198: }
0199:
0200: /** Set the whole matrix from a 1-D array
0201: * @param v 1-D array of doubles
0202: * @param columnFirst Whether to fill the column first or the row.
0203: * @throws ArrayIndexOutOfBoundsException Submatrix indices
0204: */
0205: public void setMatrix(double[] v, boolean columnFirst) {
0206: try {
0207: if (v.length != m * n)
0208: throw new IllegalArgumentException("sizes not match.");
0209: int i, j, count = 0;
0210: if (columnFirst) {
0211: for (i = 0; i < m; i++) {
0212: for (j = 0; j < n; j++) {
0213: A[i][j] = v[count];
0214: count++;
0215: }
0216: }
0217: } else {
0218: for (j = 0; j < n; j++) {
0219: for (i = 0; i < m; i++) {
0220: A[i][j] = v[count];
0221: count++;
0222: }
0223: }
0224: }
0225:
0226: } catch (ArrayIndexOutOfBoundsException e) {
0227: throw new ArrayIndexOutOfBoundsException(
0228: "Submatrix indices");
0229: }
0230: }
0231:
0232: /** Returns the maximum absolute value of all elements
0233: @return the maximum value
0234: */
0235: public double maxAbs() {
0236: double ma = Math.abs(A[0][0]);
0237: for (int j = 0; j < n; j++) {
0238: for (int i = 0; i < m; i++) {
0239: ma = Math.max(ma, Math.abs(A[i][j]));
0240: }
0241: }
0242: return ma;
0243: }
0244:
0245: /** Returns the maximum absolute value of some elements of a column,
0246: that is, the elements of A[i0:i1][j].
0247: @param i0 the index of the first element of the column
0248: @param i1 the index of the last element of the column
0249: @param j the index of the column
0250: @return the maximum value */
0251: public double maxAbs(int i0, int i1, int j) {
0252: double m = Math.abs(A[i0][j]);
0253: for (int i = i0 + 1; i <= i1; i++) {
0254: m = Math.max(m, Math.abs(A[i][j]));
0255: }
0256: return m;
0257: }
0258:
0259: /** Returns the minimum absolute value of some elements of a column,
0260: that is, the elements of A[i0:i1][j].
0261: @param i0 the index of the first element of the column
0262: @param i1 the index of the last element of the column
0263: @param column the index of the column
0264: @return the minimum value
0265: */
0266: public double minAbs(int i0, int i1, int column) {
0267: double m = Math.abs(A[i0][column]);
0268: for (int i = i0 + 1; i <= i1; i++) {
0269: m = Math.min(m, Math.abs(A[i][column]));
0270: }
0271: return m;
0272: }
0273:
0274: /** Check if the matrix is empty
0275: * @return true if the matrix is empty
0276: */
0277: public boolean isEmpty() {
0278: if (m == 0 || n == 0)
0279: return true;
0280: if (A == null)
0281: return true;
0282: return false;
0283: }
0284:
0285: /** Return a DoubleVector that stores a column of the matrix
0286: * @param j the index of the column
0287: * @return the column
0288: */
0289: public DoubleVector getColumn(int j) {
0290: DoubleVector v = new DoubleVector(m);
0291: double[] a = v.getArray();
0292: for (int i = 0; i < m; i++)
0293: a[i] = A[i][j];
0294: return v;
0295: }
0296:
0297: /** Return a DoubleVector that stores some elements of a column of the
0298: * matrix
0299: * @param i0 the index of the first element of the column
0300: * @param i1 the index of the last element of the column
0301: * @param j the index of the column
0302: * @return the DoubleVector
0303: */
0304: public DoubleVector getColumn(int i0, int i1, int j) {
0305: DoubleVector v = new DoubleVector(i1 - i0 + 1);
0306: double[] a = v.getArray();
0307: int count = 0;
0308: for (int i = i0; i <= i1; i++) {
0309: a[count] = A[i][j];
0310: count++;
0311: }
0312: return v;
0313: }
0314:
0315: /** Multiplication between a row (or part of a row) of the first matrix
0316: * and a column (or part or a column) of the second matrix.
0317: * @param i the index of the row in the first matrix
0318: * @param j0 the index of the first column in the first matrix
0319: * @param j1 the index of the last column in the first matrix
0320: * @param B the second matrix
0321: * @param l the index of the column in the second matrix
0322: * @return the result of the multiplication
0323: */
0324: public double times(int i, int j0, int j1, PaceMatrix B, int l) {
0325: double s = 0.0;
0326: for (int j = j0; j <= j1; j++) {
0327: s += A[i][j] * B.A[j][l];
0328: }
0329: return s;
0330: }
0331:
0332: /** Decimal format for converting a matrix into a string
0333: * @return the default decimal format
0334: */
0335: protected DecimalFormat[] format() {
0336: return format(0, m - 1, 0, n - 1, 7, false);
0337: }
0338:
0339: /** Decimal format for converting a matrix into a string
0340: * @param digits the number of digits
0341: * @return the decimal format
0342: */
0343: protected DecimalFormat[] format(int digits) {
0344: return format(0, m - 1, 0, n - 1, digits, false);
0345: }
0346:
0347: /** Decimal format for converting a matrix into a string
0348: * @param digits the number of digits
0349: * @param trailing
0350: * @return the decimal format
0351: */
0352: protected DecimalFormat[] format(int digits, boolean trailing) {
0353: return format(0, m - 1, 0, n - 1, digits, trailing);
0354: }
0355:
0356: /** Decimal format for converting a matrix into a string
0357: * @param i0
0358: * @param i1
0359: * @param j
0360: * @param digits the number of digits
0361: * @param trailing
0362: * @return the decimal format
0363: */
0364: protected DecimalFormat format(int i0, int i1, int j, int digits,
0365: boolean trailing) {
0366: FlexibleDecimalFormat df = new FlexibleDecimalFormat(digits,
0367: trailing);
0368: df.grouping(true);
0369: for (int i = i0; i <= i1; i++)
0370: df.update(A[i][j]);
0371: return df;
0372: }
0373:
0374: /** Decimal format for converting a matrix into a string
0375: * @param i0
0376: * @param i1
0377: * @param j0
0378: * @param j1
0379: * @param trailing
0380: * @param digits the number of digits
0381: * @return the decimal format
0382: */
0383: protected DecimalFormat[] format(int i0, int i1, int j0, int j1,
0384: int digits, boolean trailing) {
0385: DecimalFormat[] f = new DecimalFormat[j1 - j0 + 1];
0386: for (int j = j0; j <= j1; j++) {
0387: f[j] = format(i0, i1, j, digits, trailing);
0388: }
0389: return f;
0390: }
0391:
0392: /**
0393: * Converts matrix to string
0394: *
0395: * @return the matrix as string
0396: */
0397: public String toString() {
0398: return toString(5, false);
0399: }
0400:
0401: /**
0402: * Converts matrix to string
0403: *
0404: * @param digits number of digits after decimal point
0405: * @param trailing true if trailing zeros are padded
0406: * @return the matrix as string
0407: */
0408: public String toString(int digits, boolean trailing) {
0409:
0410: if (isEmpty())
0411: return "null matrix";
0412:
0413: StringBuffer text = new StringBuffer();
0414: DecimalFormat[] nf = format(digits, trailing);
0415: int numCols = 0;
0416: int count = 0;
0417: int width = 80;
0418: int lenNumber;
0419:
0420: int[] nCols = new int[n];
0421: int nk = 0;
0422: for (int j = 0; j < n; j++) {
0423: lenNumber = nf[j].format(A[0][j]).length();
0424: if (count + 1 + lenNumber > width - 1) {
0425: nCols[nk++] = numCols;
0426: count = 0;
0427: numCols = 0;
0428: }
0429: count += 1 + lenNumber;
0430: ++numCols;
0431: }
0432: nCols[nk] = numCols;
0433:
0434: nk = 0;
0435: for (int k = 0; k < n;) {
0436: for (int i = 0; i < m; i++) {
0437: for (int j = k; j < k + nCols[nk]; j++)
0438: text.append(" " + nf[j].format(A[i][j]));
0439: text.append("\n");
0440: }
0441: k += nCols[nk];
0442: ++nk;
0443: text.append("\n");
0444: }
0445:
0446: return text.toString();
0447: }
0448:
0449: /** Squared sum of a column or row in a matrix
0450: * @param j the index of the column or row
0451: * @param i0 the index of the first element
0452: * @param i1 the index of the last element
0453: * @param col if true, sum over a column; otherwise, over a row
0454: * @return the squared sum
0455: */
0456: public double sum2(int j, int i0, int i1, boolean col) {
0457: double s2 = 0;
0458: if (col) { // column
0459: for (int i = i0; i <= i1; i++)
0460: s2 += A[i][j] * A[i][j];
0461: } else {
0462: for (int i = i0; i <= i1; i++)
0463: s2 += A[j][i] * A[j][i];
0464: }
0465: return s2;
0466: }
0467:
0468: /** Squared sum of columns or rows of a matrix
0469: * @param col if true, sum over columns; otherwise, over rows
0470: * @return the squared sum
0471: */
0472: public double[] sum2(boolean col) {
0473: int l = col ? n : m;
0474: int p = col ? m : n;
0475: double[] s2 = new double[l];
0476: for (int i = 0; i < l; i++)
0477: s2[i] = sum2(i, 0, p - 1, col);
0478: return s2;
0479: }
0480:
0481: /** Constructs single Householder transformation for a column
0482: *
0483: @param j the index of the column
0484: @param k the index of the row
0485: @return d and q
0486: */
0487: public double[] h1(int j, int k) {
0488: double dq[] = new double[2];
0489: double s2 = sum2(j, k, m - 1, true);
0490: dq[0] = A[k][j] >= 0 ? -Math.sqrt(s2) : Math.sqrt(s2);
0491: A[k][j] -= dq[0];
0492: dq[1] = A[k][j] * dq[0];
0493: return dq;
0494: }
0495:
0496: /** Performs single Householder transformation on one column of a matrix
0497: *
0498: @param j the index of the column
0499: @param k the index of the row
0500: @param q q = - u'u/2; must be negative
0501: @param b the matrix to be transformed
0502: @param l the column of the matrix b
0503: */
0504: public void h2(int j, int k, double q, PaceMatrix b, int l) {
0505: double s = 0, alpha;
0506: for (int i = k; i < m; i++)
0507: s += A[i][j] * b.A[i][l];
0508: alpha = s / q;
0509: for (int i = k; i < m; i++)
0510: b.A[i][l] += alpha * A[i][j];
0511: }
0512:
0513: /** Constructs the Givens rotation
0514: * @param a
0515: * @param b
0516: * @return a double array that stores the cosine and sine values
0517: */
0518: public double[] g1(double a, double b) {
0519: double cs[] = new double[2];
0520: double r = Maths.hypot(a, b);
0521: if (r == 0.0) {
0522: cs[0] = 1;
0523: cs[1] = 0;
0524: } else {
0525: cs[0] = a / r;
0526: cs[1] = b / r;
0527: }
0528: return cs;
0529: }
0530:
0531: /** Performs the Givens rotation
0532: * @param cs a array storing the cosine and sine values
0533: * @param i0 the index of the row of the first element
0534: * @param i1 the index of the row of the second element
0535: * @param j the index of the column
0536: */
0537: public void g2(double cs[], int i0, int i1, int j) {
0538: double w = cs[0] * A[i0][j] + cs[1] * A[i1][j];
0539: A[i1][j] = -cs[1] * A[i0][j] + cs[0] * A[i1][j];
0540: A[i0][j] = w;
0541: }
0542:
0543: /** Forward ordering of columns in terms of response explanation. On
0544: * input, matrices A and b are already QR-transformed. The indices of
0545: * transformed columns are stored in the pivoting vector.
0546: *
0547: *@param b the PaceMatrix b
0548: *@param pvt the pivoting vector
0549: *@param k0 the first k0 columns (in pvt) of A are not to be changed
0550: **/
0551: public void forward(PaceMatrix b, IntVector pvt, int k0) {
0552: for (int j = k0; j < Math.min(pvt.size(), m); j++) {
0553: steplsqr(b, pvt, j, mostExplainingColumn(b, pvt, j), true);
0554: }
0555: }
0556:
0557: /** Returns the index of the column that has the largest (squared)
0558: * response, when each of columns pvt[ks:] is moved to become the
0559: * ks-th column. On input, A and b are both QR-transformed.
0560: *
0561: * @param b response
0562: * @param pvt pivoting index of A
0563: * @param ks columns pvt[ks:] of A are to be tested
0564: * @return the index of the column
0565: */
0566: public int mostExplainingColumn(PaceMatrix b, IntVector pvt, int ks) {
0567: double val;
0568: int[] p = pvt.getArray();
0569: double ma = columnResponseExplanation(b, pvt, ks, ks);
0570: int jma = ks;
0571: for (int i = ks + 1; i < pvt.size(); i++) {
0572: val = columnResponseExplanation(b, pvt, i, ks);
0573: if (val > ma) {
0574: ma = val;
0575: jma = i;
0576: }
0577: }
0578: return jma;
0579: }
0580:
0581: /** Backward ordering of columns in terms of response explanation. On
0582: * input, matrices A and b are already QR-transformed. The indices of
0583: * transformed columns are stored in the pivoting vector.
0584: *
0585: * A and b must have the same number of rows, being the (pseudo-)rank.
0586: *
0587: * @param b PaceMatrix b
0588: * @param pvt pivoting vector
0589: * @param ks number of QR-transformed columns; psuedo-rank of A
0590: * @param k0 first k0 columns in pvt[] are not to be ordered.
0591: */
0592: public void backward(PaceMatrix b, IntVector pvt, int ks, int k0) {
0593: for (int j = ks; j > k0; j--) {
0594: steplsqr(b, pvt, j, leastExplainingColumn(b, pvt, j, k0),
0595: false);
0596: }
0597: }
0598:
0599: /** Returns the index of the column that has the smallest (squared)
0600: * response, when the column is moved to become the (ks-1)-th
0601: * column. On input, A and b are both QR-transformed.
0602: *
0603: * @param b response
0604: * @param pvt pivoting index of A
0605: * @param ks psudo-rank of A
0606: * @param k0 A[][pvt[0:(k0-1)]] are excluded from the testing.
0607: * @return the index of the column
0608: */
0609: public int leastExplainingColumn(PaceMatrix b, IntVector pvt,
0610: int ks, int k0) {
0611: double val;
0612: int[] p = pvt.getArray();
0613: double mi = columnResponseExplanation(b, pvt, ks - 1, ks);
0614: int jmi = ks - 1;
0615: for (int i = k0; i < ks - 1; i++) {
0616: val = columnResponseExplanation(b, pvt, i, ks);
0617: if (val <= mi) {
0618: mi = val;
0619: jmi = i;
0620: }
0621: }
0622: return jmi;
0623: }
0624:
0625: /** Returns the squared ks-th response value if the j-th column becomes
0626: * the ks-th after orthogonal transformation. A[][pvt[ks:j]] (or
0627: * A[][pvt[j:ks]], if ks > j) and b[] are already QR-transformed
0628: * on input and will remain unchanged on output.
0629: *
0630: * More generally, it returns the inner product of the corresponding
0631: * row vector of the response PaceMatrix. (To be implemented.)
0632: *
0633: *@param b PaceMatrix b
0634: *@param pvt pivoting vector
0635: *@param j the column A[pvt[j]][] is to be moved
0636: *@param ks the target column A[pvt[ks]][]
0637: *@return the squared response value */
0638: public double columnResponseExplanation(PaceMatrix b,
0639: IntVector pvt, int j, int ks) {
0640: /* Implementation:
0641: *
0642: * If j == ks - 1, returns the squared ks-th response directly.
0643: *
0644: * If j > ks -1, returns the ks-th response after
0645: * Householder-transforming the j-th column and the response.
0646: *
0647: * If j < ks - 1, returns the ks-th response after a sequence of
0648: * Givens rotations starting from the j-th row. */
0649:
0650: int k, l;
0651: double[] xxx = new double[n];
0652: int[] p = pvt.getArray();
0653: double val;
0654:
0655: if (j == ks - 1)
0656: val = b.A[j][0];
0657: else if (j > ks - 1) {
0658: int jm = Math.min(n - 1, j);
0659: DoubleVector u = getColumn(ks, jm, p[j]);
0660: DoubleVector v = b.getColumn(ks, jm, 0);
0661: val = v.innerProduct(u) / u.norm2();
0662: } else { // ks > j
0663: for (k = j + 1; k < ks; k++)
0664: // make a copy of A[j][]
0665: xxx[k] = A[j][p[k]];
0666: val = b.A[j][0];
0667: double[] cs;
0668: for (k = j + 1; k < ks; k++) {
0669: cs = g1(xxx[k], A[k][p[k]]);
0670: for (l = k + 1; l < ks; l++)
0671: xxx[l] = -cs[1] * xxx[l] + cs[0] * A[k][p[l]];
0672: val = -cs[1] * val + cs[0] * b.A[k][0];
0673: }
0674: }
0675: return val * val; // or inner product in later implementation???
0676: }
0677:
0678: /**
0679: * QR transformation for a least squares problem<br/>
0680: * A x = b<br/>
0681: * implicitly both A and b are transformed. pvt.size() is the psuedo-rank of
0682: * A.
0683: *
0684: * @param b PaceMatrix b
0685: * @param pvt pivoting vector
0686: * @param k0 the first k0 columns of A (indexed by pvt) are pre-chosen.
0687: * (But subject to rank examination.)
0688: *
0689: * For example, the constant term may be reserved, in which
0690: * case k0 = 1.
0691: **/
0692: public void lsqr(PaceMatrix b, IntVector pvt, int k0) {
0693: final double TINY = 1e-15;
0694: int[] p = pvt.getArray();
0695: int ks = 0; // psuedo-rank
0696: for (int j = 0; j < k0; j++)
0697: // k0 pre-chosen columns
0698: if (sum2(p[j], ks, m - 1, true) > TINY) { // large diagonal element
0699: steplsqr(b, pvt, ks, j, true);
0700: ks++;
0701: } else { // collinear column
0702: pvt.shiftToEnd(j);
0703: pvt.setSize(pvt.size() - 1);
0704: k0--;
0705: j--;
0706: }
0707:
0708: // initial QR transformation
0709: for (int j = k0; j < Math.min(pvt.size(), m); j++) {
0710: if (sum2(p[j], ks, m - 1, true) > TINY) {
0711: steplsqr(b, pvt, ks, j, true);
0712: ks++;
0713: } else { // collinear column
0714: pvt.shiftToEnd(j);
0715: pvt.setSize(pvt.size() - 1);
0716: j--;
0717: }
0718: }
0719:
0720: b.m = m = ks; // reset number of rows
0721: pvt.setSize(ks);
0722: }
0723:
0724: /** QR transformation for a least squares problem <br/>
0725: * A x = b <br/>
0726: * implicitly both A and b are transformed. pvt.size() is the psuedo-rank of A.
0727: *
0728: * @param b PaceMatrix b
0729: * @param pvt pivoting vector
0730: * @param k0 the first k0 columns of A (indexed by pvt) are pre-chosen.
0731: * (But subject to rank examination.)
0732: *
0733: * For example, the constant term may be reserved, in which
0734: * case k0 = 1.
0735: **/
0736: public void lsqrSelection(PaceMatrix b, IntVector pvt, int k0) {
0737: int numObs = m; // number of instances
0738: int numXs = pvt.size();
0739:
0740: lsqr(b, pvt, k0);
0741:
0742: if (numXs > 200 || numXs > numObs) { // too many columns.
0743: forward(b, pvt, k0);
0744: }
0745: backward(b, pvt, pvt.size(), k0);
0746: }
0747:
0748: /**
0749: * Sets all diagonal elements to be positive (or nonnegative) without
0750: * changing the least squares solution
0751: * @param Y the response
0752: * @param pvt the pivoted column index
0753: */
0754: public void positiveDiagonal(PaceMatrix Y, IntVector pvt) {
0755:
0756: int[] p = pvt.getArray();
0757: for (int i = 0; i < pvt.size(); i++) {
0758: if (A[i][p[i]] < 0.0) {
0759: for (int j = i; j < pvt.size(); j++)
0760: A[i][p[j]] = -A[i][p[j]];
0761: Y.A[i][0] = -Y.A[i][0];
0762: }
0763: }
0764: }
0765:
0766: /** Stepwise least squares QR-decomposition of the problem
0767: * A x = b
0768: @param b PaceMatrix b
0769: @param pvt pivoting vector
0770: @param ks number of transformed columns
0771: @param j pvt[j], the column to adjoin or delete
0772: @param adjoin to adjoin if true; otherwise, to delete */
0773: public void steplsqr(PaceMatrix b, IntVector pvt, int ks, int j,
0774: boolean adjoin) {
0775: final int kp = pvt.size(); // number of columns under consideration
0776: int[] p = pvt.getArray();
0777:
0778: if (adjoin) { // adjoining
0779: int pj = p[j];
0780: pvt.swap(ks, j);
0781: double dq[] = h1(pj, ks);
0782: int pk;
0783: for (int k = ks + 1; k < kp; k++) {
0784: pk = p[k];
0785: h2(pj, ks, dq[1], this , pk);
0786: }
0787: h2(pj, ks, dq[1], b, 0); // for matrix. ???
0788: A[ks][pj] = dq[0];
0789: for (int k = ks + 1; k < m; k++)
0790: A[k][pj] = 0;
0791: } else { // removing
0792: int pj = p[j];
0793: for (int i = j; i < ks - 1; i++)
0794: p[i] = p[i + 1];
0795: p[ks - 1] = pj;
0796: double[] cs;
0797: for (int i = j; i < ks - 1; i++) {
0798: cs = g1(A[i][p[i]], A[i + 1][p[i]]);
0799: for (int l = i; l < kp; l++)
0800: g2(cs, i, i + 1, p[l]);
0801: for (int l = 0; l < b.n; l++)
0802: b.g2(cs, i, i + 1, l);
0803: }
0804: }
0805: }
0806:
0807: /** Solves upper-triangular equation <br/>
0808: * R x = b <br/>
0809: * On output, the solution is stored in b
0810: * @param b the response
0811: * @param pvt the pivoting vector
0812: * @param kp the number of the first columns involved
0813: */
0814: public void rsolve(PaceMatrix b, IntVector pvt, int kp) {
0815: if (kp == 0)
0816: b.m = 0;
0817: int i, j, k;
0818: int[] p = pvt.getArray();
0819: double s;
0820: double[][] ba = b.getArray();
0821: for (k = 0; k < b.n; k++) {
0822: ba[kp - 1][k] /= A[kp - 1][p[kp - 1]];
0823: for (i = kp - 2; i >= 0; i--) {
0824: s = 0;
0825: for (j = i + 1; j < kp; j++)
0826: s += A[i][p[j]] * ba[j][k];
0827: ba[i][k] -= s;
0828: ba[i][k] /= A[i][p[i]];
0829: }
0830: }
0831: b.m = kp;
0832: }
0833:
0834: /** Returns a new matrix which binds two matrices together with rows.
0835: * @param b the second matrix
0836: * @return the combined matrix
0837: */
0838: public PaceMatrix rbind(PaceMatrix b) {
0839: if (n != b.n)
0840: throw new IllegalArgumentException(
0841: "unequal numbers of rows.");
0842: PaceMatrix c = new PaceMatrix(m + b.m, n);
0843: c.setMatrix(0, m - 1, 0, n - 1, this );
0844: c.setMatrix(m, m + b.m - 1, 0, n - 1, b);
0845: return c;
0846: }
0847:
0848: /** Returns a new matrix which binds two matrices with columns.
0849: * @param b the second matrix
0850: * @return the combined matrix
0851: */
0852: public PaceMatrix cbind(PaceMatrix b) {
0853: if (m != b.m)
0854: throw new IllegalArgumentException(
0855: "unequal numbers of rows: " + m + " and " + b.m);
0856: PaceMatrix c = new PaceMatrix(m, n + b.n);
0857: c.setMatrix(0, m - 1, 0, n - 1, this );
0858: c.setMatrix(0, m - 1, n, n + b.n - 1, b);
0859: return c;
0860: }
0861:
0862: /** Solves the nonnegative linear squares problem. That is, <p>
0863: * <center> min || A x - b||, subject to x >= 0. </center> <p>
0864: *
0865: * For algorithm, refer to P161, Chapter 23 of C. L. Lawson and
0866: * R. J. Hanson (1974). "Solving Least Squares
0867: * Problems". Prentice-Hall.
0868: * @param b the response
0869: * @param pvt vector storing pivoting column indices
0870: * @return solution */
0871: public DoubleVector nnls(PaceMatrix b, IntVector pvt) {
0872: int j, t, counter = 0, jm = -1, n = pvt.size();
0873: double ma, max, alpha, wj;
0874: int[] p = pvt.getArray();
0875: DoubleVector x = new DoubleVector(n);
0876: double[] xA = x.getArray();
0877: PaceMatrix z = new PaceMatrix(n, 1);
0878: PaceMatrix bt;
0879:
0880: // step 1
0881: int kp = 0; // #variables in the positive set P
0882: while (true) { // step 2
0883: if (++counter > 3 * n) // should never happen
0884: throw new RuntimeException("Does not converge");
0885: t = -1;
0886: max = 0.0;
0887: bt = new PaceMatrix(b.transpose());
0888: for (j = kp; j <= n - 1; j++) { // W = A' (b - A x)
0889: wj = bt.times(0, kp, m - 1, this , p[j]);
0890: if (wj > max) { // step 4
0891: max = wj;
0892: t = j;
0893: }
0894: }
0895:
0896: // step 3
0897: if (t == -1)
0898: break; // optimum achieved
0899:
0900: // step 5
0901: pvt.swap(kp, t); // move variable from set Z to set P
0902: kp++;
0903: xA[kp - 1] = 0;
0904: steplsqr(b, pvt, kp - 1, kp - 1, true);
0905: // step 6
0906: ma = 0;
0907: while (ma < 1.5) {
0908: for (j = 0; j <= kp - 1; j++)
0909: z.A[j][0] = b.A[j][0];
0910: rsolve(z, pvt, kp);
0911: ma = 2;
0912: jm = -1;
0913: for (j = 0; j <= kp - 1; j++) { // step 7, 8 and 9
0914: if (z.A[j][0] <= 0.0) { // alpha always between 0 and 1
0915: alpha = xA[j] / (xA[j] - z.A[j][0]);
0916: if (alpha < ma) {
0917: ma = alpha;
0918: jm = j;
0919: }
0920: }
0921: }
0922: if (ma > 1.5)
0923: for (j = 0; j <= kp - 1; j++)
0924: xA[j] = z.A[j][0]; // step 7
0925: else {
0926: for (j = kp - 1; j >= 0; j--) { // step 10
0927: // Modified to avoid round-off error (which seemingly
0928: // can cause infinite loop).
0929: if (j == jm) { // step 11
0930: xA[j] = 0.0;
0931: steplsqr(b, pvt, kp, j, false);
0932: kp--; // move variable from set P to set Z
0933: } else
0934: xA[j] += ma * (z.A[j][0] - xA[j]);
0935: }
0936: }
0937: }
0938: }
0939: x.setSize(kp);
0940: pvt.setSize(kp);
0941: return x;
0942: }
0943:
0944: /** Solves the nonnegative least squares problem with equality
0945: * constraint. That is, <p>
0946: * <center> min ||A x - b||, subject to x >= 0 and c x = d. </center> <p>
0947: *
0948: * @param b the response
0949: * @param c coeficients of equality constraints
0950: * @param d constants of equality constraints
0951: * @param pvt vector storing pivoting column indices
0952: * @return the solution
0953: */
0954: public DoubleVector nnlse(PaceMatrix b, PaceMatrix c, PaceMatrix d,
0955: IntVector pvt) {
0956: double eps = 1e-10 * Math.max(c.maxAbs(), d.maxAbs())
0957: / Math.max(maxAbs(), b.maxAbs());
0958:
0959: PaceMatrix e = c.rbind(new PaceMatrix(times(eps)));
0960: PaceMatrix f = d.rbind(new PaceMatrix(b.times(eps)));
0961:
0962: return e.nnls(f, pvt);
0963: }
0964:
0965: /** Solves the nonnegative least squares problem with equality
0966: * constraint. That is, <p>
0967: * <center> min ||A x - b||, subject to x >= 0 and || x || = 1. </center>
0968: * <p>
0969: * @param b the response
0970: * @param pvt vector storing pivoting column indices
0971: * @return the solution
0972: */
0973: public DoubleVector nnlse1(PaceMatrix b, IntVector pvt) {
0974: PaceMatrix c = new PaceMatrix(1, n, 1);
0975: PaceMatrix d = new PaceMatrix(1, b.n, 1);
0976:
0977: return nnlse(b, c, d, pvt);
0978: }
0979:
0980: /** Generate matrix with standard-normally distributed random elements
0981: @param m Number of rows.
0982: @param n Number of colums.
0983: @return An m-by-n matrix with random elements. */
0984: public static Matrix randomNormal(int m, int n) {
0985: Random random = new Random();
0986:
0987: Matrix A = new Matrix(m, n);
0988: double[][] X = A.getArray();
0989: for (int i = 0; i < m; i++) {
0990: for (int j = 0; j < n; j++) {
0991: X[i][j] = random.nextGaussian();
0992: }
0993: }
0994: return A;
0995: }
0996:
0997: /**
0998: * for testing only
0999: *
1000: * @param args the commandline arguments - ignored
1001: */
1002: public static void main(String args[]) {
1003: System.out
1004: .println("================================================"
1005: + "===========");
1006: System.out
1007: .println("To test the pace estimators of linear model\n"
1008: + "coefficients.\n");
1009:
1010: double sd = 2; // standard deviation of the random error term
1011: int n = 200; // total number of observations
1012: double beta0 = 100; // intercept
1013: int k1 = 20; // number of coefficients of the first cluster
1014: double beta1 = 0; // coefficient value of the first cluster
1015: int k2 = 20; // number of coefficients of the second cluster
1016: double beta2 = 5; // coefficient value of the second cluster
1017: int k = 1 + k1 + k2;
1018:
1019: DoubleVector beta = new DoubleVector(1 + k1 + k2);
1020: beta.set(0, beta0);
1021: beta.set(1, k1, beta1);
1022: beta.set(k1 + 1, k1 + k2, beta2);
1023:
1024: System.out.println("The data set contains " + n
1025: + " observations plus " + (k1 + k2)
1026: + " variables.\n\nThe coefficients of the true model"
1027: + " are:\n\n" + beta);
1028:
1029: System.out
1030: .println("\nThe standard deviation of the error term is "
1031: + sd);
1032:
1033: System.out
1034: .println("==============================================="
1035: + "============");
1036:
1037: PaceMatrix X = new PaceMatrix(n, k1 + k2 + 1);
1038: X.setMatrix(0, n - 1, 0, 0, 1);
1039: X.setMatrix(0, n - 1, 1, k1 + k2, random(n, k1 + k2));
1040:
1041: PaceMatrix Y = new PaceMatrix(X.times(new PaceMatrix(beta))
1042: .plusEquals(randomNormal(n, 1).times(sd)));
1043:
1044: IntVector pvt = (IntVector) IntVector.seq(0, k1 + k2);
1045:
1046: /*System.out.println( "The OLS estimate (by jama.Matrix.solve()) is:\n\n" +
1047: (new PaceMatrix(X.solve(Y))).getColumn(0) );*/
1048:
1049: X.lsqrSelection(Y, pvt, 1);
1050: X.positiveDiagonal(Y, pvt);
1051:
1052: PaceMatrix sol = (PaceMatrix) Y.clone();
1053: X.rsolve(sol, pvt, pvt.size());
1054: DoubleVector betaHat = sol.getColumn(0).unpivoting(pvt, k);
1055: System.out
1056: .println("\nThe OLS estimate (through lsqr()) is: \n\n"
1057: + betaHat);
1058:
1059: System.out
1060: .println("\nQuadratic loss of the OLS estimate (||X b - X bHat||^2) = "
1061: + (new PaceMatrix(X.times(new PaceMatrix(beta
1062: .minus(betaHat))))).getColumn(0).sum2());
1063:
1064: System.out
1065: .println("============================================="
1066: + "==============");
1067: System.out.println(" *** Pace estimation *** \n");
1068: DoubleVector r = Y.getColumn(pvt.size(), n - 1, 0);
1069: double sde = Math.sqrt(r.sum2() / r.size());
1070:
1071: System.out.println("Estimated standard deviation = " + sde);
1072:
1073: DoubleVector aHat = Y.getColumn(0, pvt.size() - 1, 0).times(
1074: 1. / sde);
1075: System.out.println("\naHat = \n" + aHat);
1076:
1077: System.out
1078: .println("\n========= Based on chi-square mixture ============");
1079:
1080: ChisqMixture d2 = new ChisqMixture();
1081: int method = MixtureDistribution.NNMMethod;
1082: DoubleVector AHat = aHat.square();
1083: d2.fit(AHat, method);
1084: System.out
1085: .println("\nEstimated mixing distribution is:\n" + d2);
1086:
1087: DoubleVector ATilde = d2.pace2(AHat);
1088: DoubleVector aTilde = ATilde.sqrt().times(aHat.sign());
1089: PaceMatrix YTilde = new PaceMatrix((new PaceMatrix(aTilde))
1090: .times(sde));
1091: X.rsolve(YTilde, pvt, pvt.size());
1092: DoubleVector betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
1093: System.out.println("\nThe pace2 estimate of coefficients = \n"
1094: + betaTilde);
1095: System.out.println("Quadratic loss = "
1096: + (new PaceMatrix(X.times(new PaceMatrix(beta
1097: .minus(betaTilde))))).getColumn(0).sum2());
1098:
1099: ATilde = d2.pace4(AHat);
1100: aTilde = ATilde.sqrt().times(aHat.sign());
1101: YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
1102: X.rsolve(YTilde, pvt, pvt.size());
1103: betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
1104: System.out.println("\nThe pace4 estimate of coefficients = \n"
1105: + betaTilde);
1106: System.out.println("Quadratic loss = "
1107: + (new PaceMatrix(X.times(new PaceMatrix(beta
1108: .minus(betaTilde))))).getColumn(0).sum2());
1109:
1110: ATilde = d2.pace6(AHat);
1111: aTilde = ATilde.sqrt().times(aHat.sign());
1112: YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
1113: X.rsolve(YTilde, pvt, pvt.size());
1114: betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
1115: System.out.println("\nThe pace6 estimate of coefficients = \n"
1116: + betaTilde);
1117: System.out.println("Quadratic loss = "
1118: + (new PaceMatrix(X.times(new PaceMatrix(beta
1119: .minus(betaTilde))))).getColumn(0).sum2());
1120:
1121: System.out
1122: .println("\n========= Based on normal mixture ============");
1123:
1124: NormalMixture d = new NormalMixture();
1125: d.fit(aHat, method);
1126: System.out.println("\nEstimated mixing distribution is:\n" + d);
1127:
1128: aTilde = d.nestedEstimate(aHat);
1129: YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
1130: X.rsolve(YTilde, pvt, pvt.size());
1131: betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
1132: System.out.println("The nested estimate of coefficients = \n"
1133: + betaTilde);
1134: System.out.println("Quadratic loss = "
1135: + (new PaceMatrix(X.times(new PaceMatrix(beta
1136: .minus(betaTilde))))).getColumn(0).sum2());
1137:
1138: aTilde = d.subsetEstimate(aHat);
1139: YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
1140: X.rsolve(YTilde, pvt, pvt.size());
1141: betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
1142: System.out.println("\nThe subset estimate of coefficients = \n"
1143: + betaTilde);
1144: System.out.println("Quadratic loss = "
1145: + (new PaceMatrix(X.times(new PaceMatrix(beta
1146: .minus(betaTilde))))).getColumn(0).sum2());
1147:
1148: aTilde = d.empiricalBayesEstimate(aHat);
1149: YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
1150: X.rsolve(YTilde, pvt, pvt.size());
1151: betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
1152: System.out
1153: .println("\nThe empirical Bayes estimate of coefficients = \n"
1154: + betaTilde);
1155:
1156: System.out.println("Quadratic loss = "
1157: + (new PaceMatrix(X.times(new PaceMatrix(beta
1158: .minus(betaTilde))))).getColumn(0).sum2());
1159:
1160: }
1161: }
|