How to use Nonlinear Least Squares Problem without Constraints

How to use Nonlinear Least Squares Problem without Constraints

Hi together,

I'm quite new with using the MKL and also not a big mathematic freak. At the moment I try to solve a Function from the type:

ai + bi * x1 = ci * x2 / (di + x3)

x1, x2, x3 are unknown, I have initial estimations for x2 and x3.

Now the questions:
1.) Can I even use the MKL for that, espacially the strnlsp* functiopns?
2.) in the init function how can I match these function to the transfer parameters and what the... is the length of an function m ?

I found some examples, but those are not really good documented - so I really need your help.

Addition: You can not offend me my critizies my knowledge at mathematic.... ;-)

Thanks for your best effort.

30 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.

Hi,Basically speaking the problem is not clear for me. Is ai,bi,ci,di,x1,x2 and x3 scalar variables? If yes, why you can't calculate x1 asx1=(ci*x2/(di+x3)-ai)/bi?With best regards,Alexander Kalinkin

And one more question: how many equations you have in your system?With best regards,Alexander Kalinkin

one more question: "I found some examples, but those are not really good documented - so I really need your help."What examples are You talking about and how we should improve the documentation?

The following comments are based on what I think is the purpose of the calculation:

Given data vectors a, b, c, d of size n >= 3, scalar coefficients x_1, x_2 and x_3 are to be found such that the n equations

a_i + b_i x_1 = c_i x_2 / (d_i + x_3), i = 1 .. n,

are satisfied in the least squares sense.

Being curious myself about the capabilities of MKL in this task, I decided to write up a short test program, adapted from the example in the MKL documentation ("Example. dtrnlsp Usage in Fortran"), with fake data made up for a, b, c, d, and reasonable termination criteria for the iterative solution. Note that the documentation contains a typo, as flagged in the source below. (Edit 2/8/2010: the typo has been corrected).

The result is

strnlsp j74 ... PASS 3.3930975E-03 3.3930975E-03

program nlfitmkl
implicit none


external j74
real, dimension(3) ::  x = (/ 3.7, 1.9, 0.17 /)
integer :: i,j,st_cr,res,n=3,m=4,iter,iter1=10,iter2=5,rci_req
real :: eps(6) = (/0.01, 0.01, 0.01, 0.01, 0.01, 0.01 /),rs=90.0
real :: r1,r2,jac_eps=1e-3,fvec(4), fjac(4,3)
integer*8 hndl
logical :: success = .FALSE.

res = strnlsp_init(hndl,n,m,x,eps,iter1,iter2,rs)
if (res /= TR_SUCCESS) then
   print *, ' Error in strnlsp_init',res
   call mkl_free_buffers()
   stop 1

RCI_req = 0
do i=1,m
   fvec(i) = 0
   do j=1,n
      fjac(i,j) = 0
   end do
end do

do while(.not. success)
   res = strnlsp_solve(hndl,fvec,fjac,rci_req)
   if (res /= TR_SUCCESS) then
      print *, ' Error in strnlsp_solve', res
      call mkl_free_buffers()
      stop 2
   select case (rci_req)
      case (-1,-2,-3,-4,-5,-6)
      success = .true.
      case (1)
         call j74(m,n,x,fvec)
      case (2)
         res = sjacobi(j74,n,m,fjac,x,jac_eps)
         if(res /= TR_SUCCESS) then
            print *, ' Error in sjacobi'
            call mkl_free_buffers()
            stop 3
   end select
end do
res=strnlsp_get(hndl,iter,st_cr,r1,r2)   ! <---- MKL doc erroneously shows r1_r2
if (res /= TR_SUCCESS) then
   print *,' Error in strnlsp_get'
   call mkl_free_buffers
   stop 4
res = strnlsp_delete(hndl)
if (res /= TR_SUCCESS) then
   print *, ' Error in strnlsp_delete'
   call mkl_free_buffers()
   stop 5
call mkl_free_buffers()
if(r2 < 4e-3) then
   print *,' strnlsp j74 ... PASS',r1,r2
   print *,' strnlsp j74 ... FAIL',r1,r2
end program nlfitmkl

subroutine j74(m,n,x,f)
implicit none
real, dimension(4) :: &
   a=  (/ -20.394,  -21.862,  -16.822,  -21.395 /), &
   b = (/   5.137,   5.365,    3.418,     5.383 /), &
   c = (/ -1.752,   -2.554,   -2.282,    -1.88  /), &
   d = (/  2.228,    2.240,    0.868,     2.25 /)
integer, intent(in) :: m,n
integer :: i
real, dimension(n), intent(in out) :: x
real, dimension(m), intent(in out) :: f
do i=1,m
end do
end subroutine j74
  1. res=strnlsp_get(hndl,iter,st_cr,r1,r2)!<----MKLdocerroneouslyshowsr1_r2

What version are talking about?


I am not sure which version this link pertains to. I had searched the Intel site, located the documentation, and bookmarked the HTML version of the reference manual, assuming that it would be updated whenever a new release of MKL was issued.

aahh, thanks a lot!I didn't check these examples -:). We will remove the typo asap.

fyi - the typo was removed..

Thanks a lot for all this help!

Many things are now clear and understood. I translated your example to my enviroment and to C and I can observe the single steps working.

But there occures an new problem. After the second Iteration (the Jacobi Matrix is been calculated the first time in 5 step, with success) - the strnlsp_solve crash with an segmentation fault or better with the output " unknown machine code instruction".

Do you have any idear what is going wrong here?

My guess is that the MKL routines are not being called with the proper argument lists. Note that in this particular instance a function pointer is among the arguments.

If you post the C source, it may be possible to locate the bug.

[cpp]    //Optimierung
    //Here I copy my local parameters to global variables for using them in the global function
    _TRNSP_HANDLE_t handle;
    int n = 3; // Count of unknowns
    int m = m_imagePointX.size()*2;  // Count of Values
    float *x;
    x = (float*) malloc (sizeof (float)*n);
    x[0] = m_k1; // Initial values
    x[1] = m_f;
    x[2] = m_Tz;
    Ipp32f eps[6];
    eps[1] = 0.005;
    eps[2] = 0.005; 
    eps[3] = 0.005; 
    eps[4] = 0.005; 
    eps[5] = 0.005; 
    eps[6] = 0.005; 
    int iter1 = 1000;
    int iter2 = 100;
    Ipp32f rs = 90;

    int res = strnlsp_init(&handle, &n, &m, x, eps, &iter1, &iter2, &rs);
    if (res != TR_SUCCESS)
        std::cerr<<"Fehler bei Least Square Init ="< 

and here the Output on the screen:

Erfolg bei Least Square Init =1501 n3 m140
Neue Runde = 0
Erfolg bei Least Square Solve Iteration=0 RCI_Request=1
CalculateFunctionTsai - Start
CalculateFunctionTsai - End mit x0=0 x1=0.394556 x2=-0.0816883
Neue Runde = 1
Erfolg bei Least Square Solve Iteration=1 RCI_Request=2
CalculateFunctionTsai - Start
CalculateFunctionTsai - End mit x0=27 x1=0.394556 x2=-0.0816883
CalculateFunctionTsai - Start
CalculateFunctionTsai - End mit x0=-27 x1=0.394556 x2=-0.0816883
CalculateFunctionTsai - Start
CalculateFunctionTsai - End mit x0=0 x1=27.3946 x2=-0.0816883
CalculateFunctionTsai - Start
CalculateFunctionTsai - End mit x0=0 x1=-26.6054 x2=-0.0816883
CalculateFunctionTsai - Start
CalculateFunctionTsai - End mit x0=0 x1=0.394556 x2=26.9183
CalculateFunctionTsai - Start
CalculateFunctionTsai - End mit x0=0 x1=0.394556 x2=-27.0817
Erfolg bei Jacobi Iteration=1 n=3 m=140
Neue Runde = 2
Ungltiger Maschinenbefehl

The last means SIGILL

I checked this example also and have some questions about:

1.) in Line 36 (C-Code Example) MKL)INT - I think this should be MKL_INT
2.) I'm confused about the access to the elemts of x:

    for (i = 0; i < n/4; i++)


        x [4*i]     = 3.0;

        x [4*i + 1] = -1.0;

        x [4*i + 2] = 0.0;

        x [4*i + 3] = 1.0;


This will always steps through the loop once, right? Is there a deeper sense-which I do not understand?

And at the end the same again:

void extended_powell (MKL_INT *m, MKL_INT *n, double *x, double *f)


    MKL_INT i;


    for (i = 0; i < (*n)/4; i++)


       f [4*i] = x [4*i] + 10.0*x [4*i + 1];

       f [4*i + 1] = 2.2360679774998*(x [4*i + 2] - x [4*i + 3]);

       f [4*i + 2] = (x [4*i + 1] - 2.0*x [4*i + 2])*(x [4*i + 1] - 2.0*x [4*i + 2]);

       f [4*i + 3] = 3.1622776601684*(x [4*i] - x [4*i + 3])*(x [4*i] - x [4*i + 3]);




Gentlemen,First of all thank you for so deep interest in Trust Region functionality! Nevertheless Jeremy could you provide us your example fully (subroutine CalculateFunctionTsai and other) to reproduce your problem.With best regards,Alexander Kalinkin

Aboutextendet_powell function: you can set m=n=8 for example than your example will step through loop twice.With best regards,Alexander Kalinkin

Sorry that I can not share the whole code, but for a better understanding I can provide all the inputs and the code parts which are interessting for this error:

1.) The CalculateFunctionTsai Subroutine:

[cpp]void CalculateFunctionTsai(int *m, int *n, double *x, double *FValues)
	std::cerr<<"CalculateFunctionTsai - Start m="<< *m< 
2.) The Start Parameters:

[cpp]IS::Error CameraCalib::ReadCalib()

	//Kalibrierwerte ToDo: Wird bergeben
	int countSpan = 7;
	int countRow = 10;
	int countTotal = 70;
	Ipp32f xRefOffset = 100.0;
	Ipp32f yRefOffset = 100.0;
	Ipp32f xRefStepSize = 40.0;
	Ipp32f yRefStepSize = 80.0;

	Ipp32f xImg[70] ={ 42.0, 127.0, 214.0, 305.0, 399.0, 490.5, 576.0,
							   36.0, 121.5, 210.0, 304.5, 399.0, 492.0, 580.5,
							   31.5, 115.5, 207.0, 301.5, 397.5, 492.0, 580.5,
							   27.0, 114.0, 204.0, 300.0, 396.0, 492.0, 580.5,
							   22.5, 111.0, 202.5, 298.5, 394.5, 490.5, 580.5,
							   21.0, 108.0, 201.0, 295.5, 393.0, 389.0, 577.5,
							   18.0, 106.5, 199.5, 294.0, 391.5, 487.5, 576.0,
							   16.5, 105.0, 198.0, 292.5, 388.5, 484.5, 573.0,
							   16.5, 103.5, 196.5, 292.5, 387.0, 481.5, 570.0,
							   18.0, 103.5, 196.5, 291.0, 384.0, 477.0, 565.5};

	Ipp32f yImg[70] ={ 34.0,  33.0,  29.0,  31.0,  35.0,  40.5,  49.5,
							   75.0,  73.5,  73.5,  73.5,  78.0,  82.5,  90.0,
							  115.5, 117.0, 117.0, 117.5, 121.5, 126.0, 135.0,
							  162.0, 162.0, 163.5, 166.5, 165.5, 174.0, 178.5,
							  207.0, 208.5, 211.5, 214.5, 219.0, 222.0, 226.5,
							  253.5, 255.0, 259.5, 264.0, 267.0, 271.5, 273.0,
							  300.0, 303.0, 307.5, 312.0, 319.0, 318.0, 319.5,
							  346.5, 351.0, 357.0, 360.0, 364.5, 366.0, 367.5,
							  391.5, 399.0, 402.0, 408.0, 412.5, 412.5, 411.5,
							  436.5, 445.5, 451.5, 456.0, 459.0, 459.0, 456.0};


	int t=0;

	for(int y=0; y 

3.) The already calculated parameters:
m_Ty =-781.346
m_Tx =-0.772715
m_sx =1.5971
m_R1=-0.874132 m_R2=0.025665 m_R3=0.48501
m_R4=-0.407516 m_R5=-0.582054 m_R6=-0.703665
m_R7=0.264242 m_R8=-0.812745 m_R9=0.519251

I guess this is all.

There is are some informations I want to share with you:
  • I found out, that when I overwrite the results of fjac again width 0.0 after calling sjacobi(CalculateFunctionTsai,&n,&m,fjac,x,eps); (line 67)-> the error don't occures! It seems like it have to do width the content of fjac. -> I'm really confused
  • It doesn't make any difference if I use the double functions instead of the Real Functions.

Last but not least all variables called g_* are only copies of the corresponding m_* variable

There is no difference, when I limit m to n.

really difficult to reproduce your testcase from parts of program. Could you
combine really small example that I could compile and execute on my side? You
dont need to share whole code but without testcase it is really hard to find
out rootcase of your problem.With best regards,Alexander Kalinkin


I transfer all to a small C-Programm. This should help to reproduce the error:

 * test.cpp
 *  Created on: 10.02.2011
 *      Author: jeremy


	//! Innere Kalibrimain.o:(.bss+0x88): first defined here
	float m_f, m_cx, m_cy, m_k1, m_k2;

	//! Rotationsmatrix der Kamera
	float m_R1, m_R2, m_R3, m_R4, m_R5, m_R6, m_R7, m_R8, m_R9;

	float m_Tx, m_Ty, m_Tz;

	//! Sensorwerte
	float m_dxs, m_dy, m_Ncx, m_Nfx, m_sx;

	//! Vector Referenzpunkte
	std::vector m_refPointX, m_refPointY, m_refPointZ;

	//! Vector Bildpunkte
	std::vector m_imagePointX, m_imagePointY;

void CalculateFunctionTsai(int *m, int *n, double *x, double *FValues)
	std::cerr<<"CalculateFunctionTsai - Start m="<< *m< best regards Jeremy

Jeremy, may be I missed something the part of discussions,but you mentioned above that some crash happened: "the strnlsp_solve crash with an segmentation fault or better with the output " unknown machine code instruction"."I quickly checked your example - the test finished without run-time problems..

what a surprise - now I'm really confused the same code crash with the german output on the screen: "Ungltiger Maschinenbefehl" what is equal with SIGILL. Could it be that it have something to do with SIMD? I'm using a Pentium M 1500MHz, which can only provide MMX, SSE and SSE2 and SIGILL can have some reason in this...- just an idear..

We test the program on an second PC - with the same bad result. This PC has an Intel Pentium Dual T3200 CPU. So it seems to depend to other things..

need more details like these:- OS
- MKL version- is it 32 or 64 bit?- the linking line

Which Libs do you include? Could it be, that we have to search the reason in the Makefile?

This is my Makefile:

# Makefile for building: test
# Generated by qmake (2.01a) (Qt 4.7.1) on: Do. Feb 10 09:10:03 2011
# Project:
# Template: app
# Command: /usr/local/Trolltech/Qt-4.7.1/bin/qmake -o Makefile

####### Compiler, tools and options

CC            = gcc
CXX           = g++
CFLAGS        = -pipe -O2 -D_REENTRANT $(DEFINES)
INCPATH       = -I/usr/local/Trolltech/Qt-4.7.1/mkspecs/linux-g++ -I. -I/usr/local/Trolltech/Qt-4.7.1/include/QtCore -I/usr/local/Trolltech/Qt-4.7.1/include/QtGui -I/usr/local/Trolltech/Qt-4.7.1/include -I/opt/intel/composerxe-2011.2.137/mkl/include -I.
LINK          = g++
LFLAGS        = -Wl,-O1 -Wl,-rpath,/usr/local/Trolltech/Qt-4.7.1/lib
LIBS          = $(SUBLIBS)  -L/usr/local/Trolltech/Qt-4.7.1/lib -L/opt/intel/lib/ia32 -L/opt/intel/ipp/lib/ia32 -L/opt/intel/composerxe-2011.1.107/ipp/lib/ia32 -liomp5 -L/opt/intel/mkl/lib/ia32 -lmkl_sequential -lmkl_core -lmkl_def -lmkl_solver -lmkl_avx -lmkl_gf -lmkl_gnu_thread -mkl_intel_thread -lmkl_intel -lmkl_p4 -lmkl_p4m -lmkl_p4m3 -lmkl_p4p -lmkl_rt -lQtGui -lQtCore -lpthread 
AR            = ar cqs
RANLIB        = 
QMAKE         = /usr/local/Trolltech/Qt-4.7.1/bin/qmake
TAR           = tar -cf
COMPRESS      = gzip -9f
COPY          = cp -f
SED           = sed
COPY_DIR      = $(COPY) -r
STRIP         = strip
INSTALL_FILE  = install -m 644 -p
INSTALL_PROGRAM = install -m 755 -p
DEL_FILE      = rm -f
SYMLINK       = ln -f -s
DEL_DIR       = rmdir
MOVE          = mv -f
MKDIR         = mkdir -p

####### Output directory


####### Files

SOURCES       = test.cpp 
OBJECTS       = test.o
DIST          = /usr/local/Trolltech/Qt-4.7.1/mkspecs/common/g++.conf 
DESTDIR       = 
TARGET        = test

first: all
####### Implicit rules

.SUFFIXES: .o .c .cpp .cc .cxx .C

	$(CXX) -c $(CXXFLAGS) $(INCPATH) -o "$@" "$<"

	$(CXX) -c $(CXXFLAGS) $(INCPATH) -o "$@" "$<"

	$(CXX) -c $(CXXFLAGS) $(INCPATH) -o "$@" "$<"

	$(CXX) -c $(CXXFLAGS) $(INCPATH) -o "$@" "$<"

	$(CC) -c $(CFLAGS) $(INCPATH) -o "$@" "$<"

####### Build rules

all: Makefile $(TARGET)


Makefile:  /usr/local/Trolltech/Qt-4.7.1/mkspecs/linux-g++/qmake.conf /usr/local/Trolltech/Qt-4.7.1/mkspecs/common/g++.conf 
	$(QMAKE) -o Makefile
qmake:  FORCE
	@$(QMAKE) -o Makefile

	@$(CHK_DIR_EXISTS) .tmp/test1.0.0 || $(MKDIR) .tmp/test1.0.0 
	$(COPY_FILE) --parents $(SOURCES) $(DIST) .tmp/test1.0.0/ && $(COPY_FILE) --parents test.cpp .tmp/test1.0.0/ && (cd `dirname .tmp/test1.0.0` && $(TAR) test1.0.0.tar test1.0.0 && $(COMPRESS) test1.0.0.tar) && $(MOVE) `dirname .tmp/test1.0.0`/test1.0.0.tar.gz . && $(DEL_FILE) -r .tmp/test1.0.0

	-$(DEL_FILE) *~ core *.core

####### Sub-libraries

distclean: clean
	-$(DEL_FILE) Makefile

check: first

	(cd $(QTDIR)/src/tools/moc && $(MAKE))

mocclean: compiler_moc_header_clean compiler_moc_source_clean

mocables: compiler_moc_header_make_all compiler_moc_source_make_all


####### Compile

test.o: test.cpp 
	$(CXX) -c $(CXXFLAGS) $(INCPATH) -o test.o test.cpp

####### Install

install:   FORCE

uninstall:   FORCE


- Ubuntu 10.10
- 32 Bit
- what you are asking for with linking line? - All the LIBS? If yes see the provided makefile..

smth like that:-Wl,--start-group $MKLroot/libmkl_intel.a $MKLroot/libmkl_intel_thread.a $MKLroot/libmkl_core.a -Wl,--end-group -openmp -lpthread


Please try to relink accordingly the linker adviser suggestion.

I ran your code from #18 on Win7-X64 using Intel C 11.1.065 (64 bit) and the bundled MKL, compiling with the command

S:\MKL> icl /Qmkl nlfit.cpp

The program ran to completion with no error messages on a laptop with 4GB of RAM and a dual core Intel T4300.

It would help if you compiled with the -g option and ran the code under GDB, if only to determine whether the illegal instruction is generated from your code or is in the C++ runtime or the MKL library code.

Thank you very very much.....

The reason was the wrong usage of the embedded Librarries!! Thank you for this information. Everythings work now like it should....

Leave a Comment

Please sign in to add a comment. Not a member? Join today