CODEWITHSUNDEEP
c++copencvmatrixorthogonal
dwbrito
dwbrito

Reputation: 5254

How to compute the orthonormal basis of a non square (rectangular) matrix

I need to find a way of computing the orthonormal basis for range of a matrix. In matlab this function does it.

I need to do this in c/c++ and I am actually working with OpenCV

However, I haven't found anything that provides this capability in OpenCV.

I've tried working with cvSVD, but my results aren't correct.

Any clues?

Upvotes: 1

Views: 1875

Answers (4)

Orgmir
Orgmir

Reputation: 393

This is in openCV, and it works with rectangular matrix as long as m>n, according to this paper

- (CvMat *) buildOrthonormal:(CvMat *) matrix {

    NSInteger rows = matrix->rows;
    NSInteger cols = matrix->cols;

    CvMat *Q = cvCreateMat(rows, cols, kMatrixType);
    CvMat *R = cvCreateMat(cols, cols, kMatrixType);  

    for (NSInteger k = 0; k < cols; k++) {
        cvSetReal2D(R, k, k, 0.0f);

        for (NSInteger i = 0; i < rows; i++) {
            double value = cvGetReal2D(R, k, k) + cvGetReal2D(matrix, i, k) * cvGetReal2D(matrix, i, k);
            cvSetReal2D(R, k, k, value);
        }
        cvSetReal2D(R, k, k, sqrt(cvGetReal2D(R, k, k)));    

        for (NSInteger i = 0; i < rows; i++) {
            double value = cvGetReal2D(matrix, i, k) / cvGetReal2D(R, k, k);
            cvSetReal2D(Q, i, k, value);
        }

        for (NSInteger j = k + 1; j < cols; j++) {
            cvSetReal2D(R, k, j, 0.0f);
            for (NSInteger i = 0; i < rows; i++) {
                double value = cvGetReal2D(R, k, j) + cvGetReal2D(Q, i, k) * cvGetReal2D(matrix, i, j);
                cvSetReal2D(R, k, j, value);
            }

            for (NSInteger i = 0; i < rows; i++) {
                double value = cvGetReal2D(matrix, i, j) - cvGetReal2D(R, k, j) * cvGetReal2D(Q, i, k);
                cvSetReal2D(matrix, i, j, value);
            }
        }
    }
    cvReleaseMat(&R);

    return Q;
}

Upvotes: 2

user2051331
user2051331

Reputation: 77

Matlab can generate codes.Why don't you try it??? First generate then examine and finally use it,that is all

Upvotes: 0

Cloud
Cloud

Reputation: 19333

If you need an existing toolkit/library to handle this, @PureW above has provided a valid answer. If you need to implement this function yourself, you're looking for an implementation of the Gram-Schmidt algorithm.

Here is an example problem to help you verify your code:

http://www.mia.uni-saarland.de/Teaching/NAVC-SS11/sol_c8.pdf

And here is the code (please see references for full credits). PLEASE NOTE: This example assumes that you have a set of data that is scaled decently. If you have a poorly scaled matrix, you may need to consider LU-decomposition or an appropriate pivot strategy. There are useful links on this topic in the references as well.

#include <cstdlib>
#include <iostream>
#include <math.h>
using namespace std;

// example: http://www.mia.uni-saarland.de/Teaching/NAVC-SS11/sol_c8.pdf
// page 5

double a[3][3] = {
    {1.0, 2.0, 1.0},
    {0.0, 1.0, 2.0},
    {1.0, 2.0, 0.0}
};
// any column of a is a vector

double r[3][3], q[3][3];

int main(int argc, char *argv[]) {
    int k, i, j;
    for (k=0; k<3; k++){
      r[k][k]=0; // equivalent to sum = 0
      for (i=0; i<3; i++)
        r[k][k] = r[k][k] + a[i][k] * a[i][k]; // rkk = sqr(a0k) + sqr(a1k) + sqr(a2k) 
      r[k][k] = sqrt(r[k][k]);  // ||a||
      cout << endl << "R"<<k<<k<<": " << r[k][k];

      for (i=0; i<3; i++) {
          q[i][k] = a[i][k]/r[k][k];
          cout << " q"<<i<<k<<": "<<q[i][k] << " ";
      }

      for(j=k+1; j<3; j++) {
        r[k][j]=0;
        for(i=0; i<3; i++) r[k][j] += q[i][k] * a[i][j];
        cout << endl << "r"<<k<<j<<": " <<r[k][j] <<endl;

        for (i=0; i<3; i++) a[i][j] = a[i][j] - r[k][j]*q[i][k];

        for (i=0; i<3; i++) cout << "a"<<j<<": " << a[i][j]<< " ";
      }
}

    system("PAUSE");
    return EXIT_SUCCESS;
}

References:

  1. http://www.cplusplus.com/forum/general/88888/
  2. http://www.mia.uni-saarland.de/Teaching/NAVC-SS11/sol_c8.pdf
  3. http://en.wikipedia.org/wiki/Pivot_element
  4. http://math.fullerton.edu/mathews/n2003/PivotingMod.html
  5. http://www.mathworks.com/support/solutions/en/data/1-FA9A48/index.html?solution=1-FA9A48

Upvotes: 4

PureW
PureW

Reputation: 5075

You want to look into the Gnu Scientific Library which is a nice and well-tested library building on top of the BLAS-libraries. It implements a lot of different matrix operations and is usually where I would start for linear algebra stuff. Maybe one of these would suit you?

Upvotes: 1

Related Questions